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I IS.  surplementary  notes 
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"box  method" 
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prognostic  equations 
. diagnostic  equations 


•orr  4nd  by  bloc*  nuoibor) 

time  Integration  scheme 

contour  level 

file 

buffer 

space  operator 


><^TWACT  fConllm»o  on  royoroo  <<do  if  nocoooorr  idonf</r  fp  blocit  nuaibor; 

This  program  Is  an  attempt  at  a large  scale,  multi-level  general  clrculatior 
model  of  the  atmosphere.  The  original  conception  was  to  Implement  a very  flex- 
ible program.  The  major  programs  are:  INITE  (Since  a typical  problem  will  re- 
quire more  storage  than  can  be  accomodated,  this  program  parcels  the  data  and 
codes  Information  as  to  how  It  Is  partitioned.);  GEX  (A  solution  program  to  step 
the  fields  ahead  In  time.);  and  DISPLY  (This  program  displays  the  various  fields 
which  have  been  saved  from  the  GEX  run.). 
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Appendix  C:  DISPLY 
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GEH  . 

SKPXi 

SPSEP 

SPIMT 

COMTOR 

DPFT 
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Appendix  D:  FILTER 


Appendix  E:  BERROR  codes 


This  program  is  an  attempc  at  a large  scale,  multi-level  general  cir- 
culation model  of  the  atmosphere.  The  original  conception  was  to  implement 
a very  flexible  program,  not  committed  to  any  fixed  spherical  grid,  to  any 
fixed  difference  scheme,  or  any  number  of  prognostic  variables  or  form  for 
the  equations  (although,  at  this  point  in  the  evolution  of  the  program, 
the  sec  of  primitive  equations  programmed,  were  chose  found  in  Kurihara  and 
Holloway  (1967)), 

Global  is  programmed  for  Che  (GS)  grid  system,  which  accomodates  a 
wide  variety  of  spherical  grids  including  the  Kurihara  grid.  The  solution 
program  RHS  employs  the  "box  method,"  as  developed  in  Kurihara  and  Holloway. 

Since  a typical  problem  will  require  more  storage  chan  can  be  accomo- 
dated in  core,  it  is  necessary  to  parcel  the  data,  and  so,  to  code  informa- 
tion as  to  how  it  is  partitioned.  This  is  done  in  a preliminary  program 
INITE.  Then  a solution  program,  GEX,  steps  the  fields  ahead  in  time.  For 
efficient  operation  only  chat  data,  which  is  necessary  for  updating  one 
latitude  band  with  no  computational  redundancy,  is  accomodated  in  core. 
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(11)  Kurlhara,  Yoshlo  and  Holloway,  J.  Lelch:  "Numerical  Incegraclon 
of  a nine  level  Global  Prlmlclve  Equations  Model  Formulated  by 
the  Box  Method,"  Monthly  Weather  Review:  Vol.  95,  No.  8,  August 
1967. 

(Ill)  Klein,  R.  and  Mettauer,  J.,  "The  Numerical  Solution  of  Partial 
Differential  Equations,"  January  1969. 

(Iv)  Klein,  Mettauer,  Magllone,  Spiegel:  "Certain  Finite  Difference 
Methods  for  the  Solution  of  Large  Scale  Circulation  Problems," 
Final  Report,  January  1,  1972. 

(v)  Klein  et  al. : "Certain  Finite  Difference  Methods  for  the  Solu- 
tion of  Large  Scale  Circulation  Problems  II,"  November  8,  1973. 

(vl)  Ahlberg,  Nllson,  Wash:  "The  Theory  of  Splines  and  Their  Appli- 
cations," Academic  Press,  1967. 

We  shall  hereafter  refer  to  these  papers  as 


(i) 

K paper 

(li) 

KH  paper 

(ill) 

KM  paper 

(iv) 

KMI  paper 

(V) 

KMII  paper 

(vl) 

ANW  paper. 

III.  Description  of  Model 

For  the  sake  of  completeness,  the  model  equations  below  are  taken  from 
the  KH  paper. 

The  equations  of  motion  In  X,  9,  o coordinates  are: 

- -DjCu)  + ^ u)  P.v 


J 


2 


1 


1 


f 


3(P^v) 

3c 


^ ^ s ( - tan6  „ 
-D3(v)  - [t  + uj  ?^u 


il*!  . A ^ p 

a39  3a  { * a39J  ' 


where  f is  the  Corolis  parameter;  R,  the  gas  constant  for  dry  air;  T,  tem- 
perature; t,  geopotential  of  a sigma  surface;  and  Fg,  frictional  forces  in 
the  zonal  and  meridional  directions  respectively;  and  the  three-dimen- 
sional divergence  operator  defined  in  KH.  The  first  law  of  thermodynamics 


3t 


-D.(T)  + — — + 

3 *^9  ^p  ^ 


where  Cp  is  the  specific  heat  capacity  at  constant  pressure;  u vertical  P- 
d P * 

velocity  q added  heat  per  unit  mass;  and  F~,  the  affect  of  thermal  dif- 
d t 

fusion.  The  continuity  equation  takes  the  form: 

I7*  - -f  D,(l)  ia 

where  D2  is  the  two  dimensional  divergence  operator  defined  in  KH. 

The  vertical  a and  P-velocities  are  obtained  by  the  diagnostic  relations: 


D2(1)  da  - D2(1)  do 

0 h 


“ ■ “ Hr  + » _ 


Finally,  the  hydrostatic  relation  is: 


5 - RT  - 


} 
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GLOBAL  Integrates,  in  time,  finite  difference  versions  of  these  equa- 
tions which  are  defined  on  a generalized  grid  (GS)^.  In  deriving  the  finite 
difference  forms  of  the  original  differential  equations,  it  is  desirable  to 
preserve  conservation  properties.  This  can  be  accomplished  by  employing  the 
"box  method"  as  an  overall  strategy.  That  is,  if 

||-  -V(^)  + F 

where  x is  a scalar  quantity;  v,  a three  dimensional  velocity  vector;  7*,  a 
three  dimensional  divergence  operator;  and  F a source  term  of  the  quantity 
x;  then  integrating  around  a volume  element  or  a box,  and  employing  the 
Gauss  Theorem  gives: 

i^ere  x is  a mean  value  of  x for  the  box  with  volume  V and  surface  S;  and 
Is  the  outward  component  of  velocity  on  the  box  surface. 

Then,  the  finite  difference  forms  of  the  prognostic  equations  become: 
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where  the  o subscript  indicates  an  average  value  of  the  variable;  taken  :c 
be  the  value  of  the  center  of  the  box,  and  where 
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and  the  space  operators  D,  H,  G;^,  Gg , L^,  Lg  are  as  defined  in  KH. 
The  diagnostic  relations  are: 

1 k k -• 

k-1  ^ k-1  ^ ^ 
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which  leads  to  the  expression  for  0 at  the  top  and  bottoia  of  a box: 


(7) 


♦lt-%  V + CTk.^) 


The  direction  of  calculation  for  each  time  step  is: 

1.  Determine  ^ recursively  at  all  levels  from  the  diagnostic  rela~ 

tlon  (7)  using  complete  knowledge  of  ~ 

2.  Calculate  and  Gg(P^). 

3.  Determine  uT  at  all  levels  from  (5)  where  the  boundary  conditions 
are  fixed  3^  ■ ^k+^j  ” 

4.  Calculate  the  tendency  from  (4). 

5.  Calculate  u>  at  all  levels  from  (6). 

6.  Calculate  the  space  operators  at  all  values  of  the  arguments 
that  appear  In  (1)  thru  (3). 

7.  Calculate  the  tendencies  of  the  3 dimensional  prognostic  varl- 
bles  from  (I)  thru  (3). 

8.  Step  the  fields  ahead  using  one  of  the  time  Integration  schemes. 

9.  Repeat  the  procedure  the  desired  number  of  time  steps. 

IV.  History 

Check  out  on  GLOBAL  began  with  a series  of  baratroplc  runs.  The  first 
significant  result  was  a leapfrog  run  on  a coarse  16  x 32  grid  and  wave  num- 
ber one;  which  went  out  to  about  60  days  and  then  blew  up.  When  the  same 
run  was  made  with  wave  number  four  the  blow  up  occurred  at  14  days  (At  was 
3608  In  this  rtin).  The  results  were  muddy  at  7 days,  and,  in  fact,  there  was 
a notlclable  tilting  in  the  v field  after  one  day.  It  seemed  that  trouble 
was  coming  down  from  the  poles,  so  it  was  decided  to  refine  the  grid.  A 
32  X 64  leapfrog,  wave  number  4 run,  with  At  ■ 75a,  ran  out  to  14  days  with- 
out trouble.  Total  energy  changed  only  in  the  9^^  decimal  place. 
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Another  approach  to  the  problem  was  to  mix  Adams  Bashford  steps  in  with 
the  leapfrog  steps.  In  the  hopes  of  stabilizing  the  run.  If  the  mixing  ratio 
was  2/100,  the  run  went  to  14  days.  However,  when  the  ratio  was  decreased 
to  1/500,  the  run  blew  up  after  3327  time  steps  (3360  ts.  = 14  days). 

To  remedy  the  disruptive  Influence  of  small  scale  effects  near  the  poles, 
various  filters  were  tried.  First,  a program  was  written  to  do  a spectral 
filter  on  the  prognostic  fields  after  every  N time  steps  (H  specified  at 
object  time).  The  program  allowed  the  user  to  specify  which  latitude  bands 
he  wished  filtered,  and  which  wave  numbers  he  wished  to  eliminate.  When 
this  filter  was  used  to  eliminate  the  2 grid  Interval  wave  from  the  height 
field,  and  from  the  winds  u,  v,  the  plots  showed  sharp  division  at  the  tran- 
sition between  unflltered  bands  (midlatitudes),  and  filtered  bands  (near  the 
poles).  To  avoid  the  disadvantage  of  such  a sharp  truncation  of  wave  number, 
a different  filtering  technique  was  applied  both  east-west  and  north-south 
on  a 16  X 32  grid.  The  filter  described  In  section  V was  applied  with  p ~ 3 
to  the  coarse  grid  after  every  10  time  steps.  The  run  resulted  In  fields 
with  little  amplitude,  and  left  unchanged  the  pronounced  tilting  apparent 
In  the  unflltered  run.  It  appears  that  most  of  the  smoothing  Is  due  to 
the  difference  scheme,  and  that  this  fact  was  disguised  In  the  unflltered 
run  by  the  blow  up. 

Another  approach,  that  was  partially  developed,  was  to  replace  the 
fields  near  the  poles  by  weighted  least  squares  polynomials,  tying  the  poly- 
nomials not  only  to  grid  values,  but  to  numerical  derivatives  near  the  poles. 
Replacing  finite  difference  forms  for  space  derivatives  near  the  poles,  with 
the  polynomials'  derivatives,  and  transforming  back  to  the  X,  6 coordinates. 

It  would  be  possible  to  prognosticate  avoiding  the  singularities  that  occur 
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in  the  finite  difference  forms  when  using  either  a polar  box  or  a latitude- 
longitude  grid. 

Fi(oJective  space  Is  defined  by: 

X ■■  x(X,6)  ■ cos6  cosX 

y ■ y(X,0)  ■ CO80  slnX  ’ 

then 

P(X,0) — --PCx.y) 


and  If  we  define: 
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Ic)  D(l)  - la  [p(l)  - P*  D(P*-10} 
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The  test  agreed  with  these  results.  Next,  the  model  was  run  with  u 
maintained  as  zero,  and  It  was  found: 

2a)  Same  as  lb). 

2b)  Same  as  Ic). 

2c)  ^ )(  0 as  expected. 

Finally,  when  the  model  was  run  with  IS  maintained  as  zero.  It  was  found: 

3«)  u - -o  P,  D(P,'*) 

3b)  ^-0 
8a 

3c)  -1^  Is  different  than  that  found  In  2c). 

9C 

V.  Filter 

The  filter  applied  east-west  Is 

> (-UP  f,  . 

where  f^  Is  the  current  variable  field,  and 
• 'i-i  - * 'i+i 
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The  filter  applied  north-south  Is: 


- [”l  + (-1)*"  gj  5 F**®  (F®”f^) 

where  gj  Is  the  field  coming  out  of  the  east  west  filter,  and:, 

"•*J  ■ <*J-1  - 'j>  - ‘fj  - Vl’  “•*1'**1 

The  above  filter  was  applied  with  p * 3,  both  east-west  and  north- 
south  on  a 16  X 32  grid.  does  not  eliminate  the  effect  of  the  2-grld 

Interval  wave.  In  fact.  If  v Is  the  two  grid  Interval  wave: 

F**®  V ■ (1  - cos^a)  V 
where  a ■ 66/2. 

Thus  F^®  V •.  .019v.  The  Legendre  polynomials  and  [sln^  are  also 
NS 

eigenvectors  of  F ; In  fact: 

- CO 

. I 1 - x‘J  51^83 

where  X ■ -2(cosa  - cos3a) 

Hethod 

After  the  completion  of  a regular  GEX  time  step,  NTS  Is  tested  to  de- 
termine If  the  current  time  step  Is  to  be  filtered.  If  so,  all  the  lOT  data 
Is  brotight  band  by  band  back  Into  BLK  (each  record  holding  384  words).  After 
this  data  Is  filtered,  it  replaces  the  original  data  on  lOT,  then.  If  INI 
data  (NTS-1)  was  not  filtered  previously.  It  Is  filtered  now. 

For  a detailed  description  of  the  filter  see  the  Appendix  D. 
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GLOBAL  Is  organized  Into  update  decks,  with  unique  deck  names  but  re-  i 

I 

dundant  program  names  within  the  decks.  This  arrangement  allows  the  modular  \ 

structure;  for  example  a program  TIS  could  be  called  from  GEX,  after  being 

compiled  within  the  deck  GTISLF.  In  this  way  a leapfrog  scheme  Is  performed 

with  the  same  call  as  would  occur  If,  Instead,  the  deck  GTISAB  were  compiled.  i 

In  that  case  the  call  to  TIS  would  execute  Adams  Bashforth  scheme.  I 

I 

The  complete  GLOBAL  program  Is  presently  In  two  card  boxes,  constituting  t 

an  update  library.  The  structure  of  the  program  organized  from  this  library 

( 

Is  determined  by  the  update  control  card,  the  update  ^COMPILE  directive  and  j 

k 

update  Insertion  or  deletion  directives.  The  box  consists  of:  I 

1 

j Comdecks  which  dimension  variables  within  programs,  and  estab- 

; llsh  common  blocks  and  blank  common  by  a descending  hierarchy  of 

1 calls  at  compile  time. 

( 

2nd  Qecks  containing  program  modules  which  may  be  organized  Into  any 
I one  of  the  three  major  types  of  runs. 

1 

I The  actual  card  file  looks  like: 

!■ 

C 

t 


i- 

t 


DISPLAY  files  saved  from 
GEX.  run 

PLOT 

REPORT 

SPECTRAL 

DUMP 


Module  options  for 
flexible  runs 


ESTABLISH  GRID 

CREATE  INITIAL  DATA  FILE 


IHITIATE  A RUN 
PROGNOSTICATE 
SAVE  DISPLAY  PILES 


DIMENSION  ARRAYS 
ESTABLISH  NAMED  COMMON 
ESTABLISH  BLANK  COMMON 


VII.  Hacro  Flowchart 


dlTE:  INITE  is  a self  contained  program  to  generate  the  data  file  and  the 
grid  conditions  for  a run. 

GEX:  Is  the  solution  program,  which  accepts  Initial  data  and  a defining 
grid  from  INITE,  and  steps  the  fields  forward  In  time.  GEX  also 
saves  user  specified  time  steps  for  later  display.  Optionally,  GEX 
may  accept  two  contiguous  time  steps  from  an  Interrupted  run,  and 
continue  that  riin.  GEX  stands  for  Global  Experiment. 

DISPLY:  DISFLY  Is  a self  contained  program  which  will  aitalyze  and  report 
the  data  generated  by  the  solution  program  GEX. 

In  the  following  we  shall  describe  all  three  modules.  However,  since 
many  of  the  variables  appear  In  two  or  all  three,  we  shall  describe  a 
variable  only  once,  the  first  time  It  appears,  to  avoid  redundancy.  Often 
this  will  result  In  describing  variables  Chat  appear,  but  are  not  used  in 
a subroutine,  and  are  only  used  later. 
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APPENDIX  A:  INITE 


Characteristics ; 

Core : 64000^ 

Tlae:  20  seconds  should  be  enough  > 

• 1 

Options:  May  choose  display  of  initial  data  field  and  grid  on  | 

community  output  file  ^ l 

Permanent  Files:  Initial  data  file 
Grid  file 


Organization  of  Core  During  a Run 


Address 

Length 

Named  Common 

/BDF/ 

111 

14600 

/UNITS/ 

14711 

15 

/AC/ 

14726 

2260 

Programs 

INITE 

17206 

6234 

REPORT 

BERROR 

OBAL 

C(»1P 

SLAY 

RINT 

PRTR 

PRTI 

PRTR2  30531  174 

CBDFN 


SYSTQI 


Blank  Common  / / 47552  12000 

Common  Blocks 

/AC/  This  block  is  used  to  store  the  variables  which  define  a 3-dimensional 
grid,  and  also,  variables  describing  the  data  organization  of  the  ISTAR  1 file. 
/UNITS/  Holds  names  of  files  that  are  manipulated  during  complete  run  of 
program. 
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INPUT  OUTPUT  CHAPT 


ISTARl 

file 


SAMPLE  DECK 


MZaZIZ— 

^^TCMA  VAT.ITF.!; 


CONST 


J^DENTIFIER  C^l  jU 


^^>^NAME,CM,T20. 


Values  of  Physical  Constants 
and  Sigma  Levels 
Parameters  for  this  run 
Identify  this  run 


Update  directive  to  compile 
specified  decks 


Update  directives — insertions 
or  deletions 


Catalog  two  files  above 


Execute 


Compile  update  file 
Organize  update  file 
Reserve  space  for  two  permanent 
files  above 
Attach  GLOBAL  library 


/BUF/  This  block  Is  used  to  hold  the  initial  data  created  in  OBAL  before 

, it  is  buffered  out  to  the  ISTARl  file.  For  the  original  goal  of  25 

levels  and  32  cells,  6528  locations  was  to  hold  exactly  two  data  bands 

of  4 three  dimensional  variables 
—levels 

2 X {(2  + 4 X 25)  X 32}  - 6528 
* t-  cells 

I diagnostic  prognostic 

variables  variables 

Files 

Ill?Tn(5) — Community  Input  File 
t OUTPUT (6) — Community  Output  File 

I 

INIT(7) — Private  Binary  Output  File  of  Grid  Information  (AC  block) 

ISTARl (8)— Initial  Data  File  Created  in  OBAL 

IDES (9) — Scratch  File — used  to  expand  data  file  if  IPT  ■ 1. 

(10) — Not  now  used 
Description  of  Variables 

> A Mesn  radius  of  the  earth 

ADT  Accumulated  delta  t.  Total  simulated  time. 

ARE  Wave  number  (R)  for  a baratropic  atmosphere  with  a free  surface. 

AVRE  Angular  rotation  of  the  earth. 

BAND  Latitude  band  boundary  array,  BAND(J)  * 9j_}j(KMI) 

BANDJ  BANDJ  ■ BAND(J)  for  current  band. 

I BCTHJ  cos(8j_ij) 

BLR  Buffer  area  for  creation  of  compressed  and  expanded  files. 

^ . CBLR  Equlvalenced  to  AC  common.  Used  for  I/O  of  (AC  block) 

[ I 

I CP  Specific  heat  at  constant  pressure  Cp. 

i 

f 


CTHJ  Array  holding  cos(6j) 

DEV  Array  holding  deviations  of  cell  (ItJ)  from  the  zero  meridian  in 
fractions  of  AXj.  If  y*  is  the  angle  of  deviation,  then 
DEV(J)  - Oj  - Yj/AXj  (KSI). 

DL  Array  of  zonal  angle  increments  DL(J)  * AXj. 

DLL  Meridional  angle  increment  DLL  ■ ®J-Jj  " ®j+is* 

DT  Simulated  time  increment  At.  Value  is  read  in  as  parameter  on  entering 
GEX. 


DTD  2At 
DTH  At/2 

DTHJ  - 9j+;j)  -Ay2  5 a 

DTI  Duration  of  the  time  integration  in  simulated  time. 

DUMS  Array  of  AW  ■ Wjj  - Wg,  where  Wj^  ■ Z and  Wg  ■ .1  ; then  DWNS(J)  ■ AWj. 

N S 

Used  in  calculation  of  Gg  operator  (KH) . 

DX  Array  of  cell  lengths  at  band  boundaries.  DX(J)  ■ a cos(0j_;j)  AXj. 

This  is  a temporary  array  and  li.  not  stored  on  the  INIT  tape  of  grid 
information. 

DXG  Array  of  cell  lengths  at  grid  points.  DXG(J)  • a coa(9j)  AXj.  This 
is  a temporary  array  and  is  not  stored  on  the  INIT  tape  of  grid 
information. 

DY  Array  of  widths  of  cells  on  band  J.  DY(J)  - aAOj.  This  is  a temporary 
array  and  is  not  stored  on  the  INIT  tape  of  grid  information. 

FIT  Fixed  interval  of  time.  Simulated  time  increment  At  in  GEX  run. 


FMB 

G^ 

HFT 


Array  of  floated  values  of  bj . 
Gravitational  constant  g. 

- 1 fff 

Array  of  values  of  f ■ — j f dV 


2Gsln0j^jCOsC^®iJ/2}  (KH). 


AV 


HMT  Array  of  metric  term  coefficients.  .®  * 1 I « tan044/a  (KH). 


AV 


'ij' 
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IDC  60  character  Identifier  for  INITE  run. 

IDES  File  name  for  display  file  catalogued  as  a permanent  file  In  GEX  run. 
Used  only  as  a scratch  file  In  an  INITE  run. 

IDESF  Frequency  of  display  file  generation  In  a GEX  run. 

JMAX  Zonal  Index  of  maximum  advectlng  velocity. 

INIT  Private  binary  output  file  of  grid  Information  saved  In  AC  block. 
Catalogued  as  a permanent  file  In  an  ItilTE  run. 

INPUT  Community  Input  file. 

INI  Expanded  active  data  file  for  time  step  N. 

IN2  Expanded  active  data  file  for  time  step  N - 1 

lOSZ  Size  of  I/O  buffer  BLK.  Presently  holds  2 full  bands  of  data  or 
6528  words. 

lOT  Expanded  active  data  Jlle  for  time  step  N + 1. 

IPRT  Community  output  file. 

IPT  Print  option  for  INITE  run. 

IPT  “0  No  point 

IPT  ■ 1 Print  all  grid  Infor^tlon  and  data 
ISC  Scratch  file.  Not  used. 


ISER  Equlvalenced  to  BLK(l).  Used  to  pass  serial  number  In  Integer  format 
off  label. 

ISS  Code  for  abnormal  stop  If  P^^  < 0. 

ISS  ■ 0 No  stop 

ISS  ■ 1 Stop 

ISTARl  Initial  data  file  for  INITE  run.  Initial  data  or  restart  data  file 
for  subsequent  GEX  run. 

ISTAR2  Restart  data  file  for  GEX  run. 

ITS  Equlvalenced  to  BLK(2).  Used  to  pass  time  step  number  In  Integer 
format  off  label. 

II  INIT  file  name.  Avoids  technicality  that  can't  pass  file  name  In 
common  to  unit  function. 

J Latitude  band  Index. 
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JB  Number  of  data  bands  In  grid.  Set  as  one  of  the  parameters  read  In 
INITE  run. 

JBl  JB  + 1.  Actual  Index  of  last  data  band,  since  first  data  band  at 
- Nocth.  jPole  is  a dummy. 

JBIZ  JBl.  Avoids  use  of  JB  in  common  as  control  on  Do  loop. 

JB2  JB  + 2.  Index  of  first  dummy  band  at  South  Pole. 

JB3  JB  + 3.  Index  of  second  dummy  band  at  South  Pole. 

JB4  JB  + 4.  Index  of  third  dummy  band  at  South  Pole. 

JBS  JB  + 5.  Index  of  last  dummy  band  at  South  Pole. 

JLI  Array  of  indices  of  the  last  bands  fit  into  the  blocks  of  the  compressed 
data  file. 

JMAX  Latitude  index  of  the  maximum  advecting  velocity. 

JMIN  Latitude  index  of  minimum  space  increment. 

JX  Maximum  number  of  data  bands  that  can  be  accomodated. 

K Vertical  level  index. 

KMAX  Vertical  index  of  maximum  advecting  velocity. 

KMIN  Level  index  of  minimum  space  Increment. 

KV  Number  of  vertical  levels. 

KVM  KV  - 1. 

KV2  KV  + 1.  For  KV  levels,  there  are  KVl  boundaries. 

KX  Maximum  number  of  levels  that  can  be  accomodated.  This  is  limited 

by  the  second  dimension  of  the  P arrays.  Presently,  P(4,  3,  32),  so 
KX  - 3 - 2 - 1. 

K1  Storage  location  in  P arrays  for  diagnostic  variables  and  grid  infor- 
mation. 

K2  Storage  location  in  P arrays  for  diagnostic  variables  and  grid  infor- 
mation. 

L Index  for  arrays  in  INITE  run. 

LO  Smallest  time  step  number  on  a display  file.  Used  to  seek  file  in 

display  run. 
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Largest  tine  step  number  on  a display  file.  Used  to  seek  file  In 
display  run. 

Number  of  cells  on  each  band  (fixed)  , or.  If  LLR  • 0 code  to  Indicate 
that  a variable  number  of  cells  on  each  band  (bj)  must  be  read  In. 

Accusnilated  number  of  words  of  data,  up  to  maximum  of  I/O  size.  Used 
when  partitioning  data  Into  blocks  for  I/O. 

’ J2 

Arrary  of  number  of  data  words  In  an  I/O  block.  LSZ(L)  ■ £ -LTI(J), 

J1 

where  J1  - JLI(L  - 1)  and  J2  - JLI(L). 

1 

Array  of  number  of  data  words  In  a band.  LTI(J)  > bJ  x {KV  x MP  + 2): 
that  is,  KV  levels  of  MP  3-dlmenslonal  variables,  plus  2 two  dimensional 
variables,  gives  number  of  words  required  for  each  column;  and  bj  such 
columns  around  the  band. 

Number  of  3-dlmensional  prognostic  variables.  Set  as  a read  In  para- 
meter In  INITE. 

Maximum  time  step  desired  in  a GEX  run. 

Array  of  number  of  cells  on  each  band  (bj). 

Maximum  value  In  array  NB  above. 

Number  of  records  the  compressed  data  has  been  partitioned  Into. 

Running  count  of  the  number  of  display  files  generated  in  a GEX  run. 

Number  of  data  words  required  In  a vertical  column  (KV  x MP  -t-  2). 

First  word  on  data  file  label.  A number  used  to  uniquely  Identify 
INITE  creation  run  data. 

Time  step  currently  being  processed  during  a GEX  cun. 

Gas  constant.  Don't  confuse  with  wave  number  (ARE). 

Real  time  variable. 

Real  time  limit.  Used  to  Initiate  stop  processing  In  a GEX  run.  If 
time  limit  Is  exceeded. 

Array  of  sigma  values  SIG(k)  ■ 

SIG(l)  * top  of  atmosphere 

SIG(KVl)  - bottom  of  atmosphere 

Array  of  Ao.  SIGl(k)  • ” ®k-Ji*  finite  difference  forms. 


J 


WNS  Array  of  north-south  weights. 

ZB  Array  uT.  Value  Is  zero  at  top  and  bottom  of  the  atmosphere. 

Note;  Some  confusion  may  arise  as  R Is  used  both  for  the  gas  constant,  and 
In  the  usual  notation  for  wave  number  (ARE) . 


Generate 

3rD 

GrU 


Create  A 

Latitude .Longitude 
Grid  With  No 
Band  Shlftln 


Read 
Parameters 
specifying  grid 
NB,  BAND,  DEV 


Read  the  KV  levels 
of  the  normalized 
pressure  coordinate 
SIG 


Print  # cells  on  each 
band  (bj)  Boundary 
locations  (9j  * 'j) 
Deviations  (y j ) 

Sigma  Values  (oj^) 


FILL  IN  SIGMA  TABLES 
SIGI..A0  SIG4. .RAo 

SIG2..%Ao  SIG5..  R 5 
SIG3..5  Cpd 


Fill  In  latitude  tables 
Table  of  6j,  cos^6J). 

» - A0j/2,  k,  *1  and 
weights. 


Change  band  boundaries 
to  radian  measure 
Change  deviations  to 
fractions  of  zonal 


Fill  In 

Dummy 

Slots 
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Helpful  Comments 

DLL  - 180. /FLOAT (JB)  $ BAND(l)  - 90.  $ BAND(2)  - 90.  (Ginit.61) 

There  Is  a single  band  at  the  north  pole  corresponding  to  the  Index 
J ■ 1.  Thus  the  north  pole  Is  found  at  0^  ■ 90,  where  BAND(l)  - 0^  and 
BA1ID(2)  - 03^2* 


®JW-3/2  “ ®JW-5/2  • " ®JB4-9/2  " 

JBl  Is  the  Index  of  the  last  data  band. 

JB2»  JB3,  JB4,  JB5  correspond  to  dummy  bands  at  the  south  pole,  each 
having  one  dummy  cell.  (IIB(I)  “ NB(2)  ■ NB(3)  ■ 11B(4)  • NB(5)  ^ 1) 
Actually  ire  should  be  dimensioned  69  to  accomodate  64  data  bands. 
teed  (INPDT.  102)  fSIG(k).  k - 1.  kvU (GINIT,89I 

cc.::..IIote  that  Is  Indexed  from  the  top  to  the  bottom  of  the  atmosphere, 

and,  that  SIG(k)  corresponds  to  of  the  KH  paper  and  is  the  upper  and 

lower  boundary  of  a cell. 


°1 


®l/2 

®3/2 


SIG(l)  top  atmosphere 


1 


®kv-% 


®kv4Js  siG(kvl)  bottom  of 
atmosphere 


CALL  pm  flPRT.  ARRAY.  JJ) 

Prints  a 1-dlmenslonal  Integer  array  on  the  community  output  file.  Prints 


JJ  elements  In  rows  numbering  each  row. 

Pirea)  ■ NB(J) CGINIT.  103) 

Floats  bj.  This  array  saved  on  INIT  file,  along  with  bj. 

DO  5 J ■ 2.  JBIZ  (GINIT.122) 

J • 2 thru  JBl ‘correspond  to  the  data  bands.  See  variable  description 
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Compute  number  of  Data 
blocks  (NBLK)  required  to 
haindle  data  with  I/O  size 
specified,  along  with  the 
data  size  of  these  blocks 
LSZ (L) , and  the  index  of 


Buffer  out 
INIT  file 


i 


1 

1 


Helpful  Comments 

C ifagrTbglc  Section  (GINIT.158) 

I 7 

; Here  are  created  the  descriptive  arrays  LTI,  JLI,  LSZ,  and  the  parameter 
NBLK. ''These  completely  describe  the  organization  of  the  data  file  ISTARl, 
and  all  subsequent  data  files  from  GEX  runs.  These  data  files  are  organized 
In  compressed  format — that  Is  only  the  3-D  prognostic  variables,  and  the 
variables  P*  and  are  carried;  other  variables  are  not  carried  since  they 
can  be  recreated  from  the  IHIT  file  at  the  beginning  of  each  GEX  run  (that 

t 

1®  5) • 

Organization  of  the  Data  File 

Each  data  file  consists  of  a 3 work  record  uniquely  describing  the  file; 
then  the  data  In  compressed  blocks,  each  of  size  < lOSZ.  Each  block  contains 
one  or  more  bands  of  data,  reported  column  by  column  from  the  zero^^  meridian. 

For  example,  an  Initial  data  file  (ISTARl)  with  3 dimensional  variables 
u,  V,  T would  be  organized  as  follows: 
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A 

EOF 


where 

^IJ*  '*IJK 

♦iJ*  VIJK 

I * zonal  Index  of  cell 
J * meridional  Index  of  band 
K * vertical  Index  of  O level 

1 
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r'l  ■ ' 

I 

I Helpful  Comments 

I CALL  PRTR2(IPRT.  WNS,  JB5.  2.  JX.  2) (GINIT.208) 

' Prints  a two-dimensional  real  array,  numbering  the  rows  and  columns. 

I The  parameters  In  the  third  and  fourth  positions  give  the  number  of  rows 

I and  columns,  respectively,  to  be  printed.  The  fifth  and  sixth  parameters 

f 

\ give  the  overall  dimensions  of  the  array.  In  this  case  as  dimensioned  In 

INITE  WNS (64,  2).  Notice  that  JX  has  been  set  to  64  in  the  INITE  data 


statement. 

SUBROUTINE  OBAL 

This  version  of  OBAL  was  written  specifically  for  the  baratroplc  case 
with  a free  surface.  In  particular  (1959  Phillips): 

a* 

If  c = cos6  and  u*1c«  1/100  and  R taken  as  the  wave  number 
|2 


then  I#  ■ #0 


afl 

10 


I 10. 5c^  + kc^^  QR+Dc^  + (2R^-R-2)  - 2R^c'^3 


22c^ 


(R+l)(R+2) 


[2r^+2R+2)  - (R+l)^c^n  cosRA 


+ itc^^  Qr+1)c2  - (R+2)3coa2RA} 


where  ® is  the  height  of  the  free  surface — notice  that  9 displaces 

P*- 

Hq  ■ 8 km  (top  of  atmosphere)  g^  ■ 9.80616m/s^ 
a ■ 6.371229  x 10^  m (mean  radius  of  the  earth) 
so  that 

- g H - 7.844928m^/s^ 

,-5 


0 ■ 7.292116  X 10  ^rad/s  (Coriolis  parameter) 


then 


1 

a ^ 


1 


_____  iii. 
a cosO  3\ 


1 


I 

i 
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where 


10  I 


sln6  + c slnScosRX 


1 


so  that 


anf 
“ " lOL 


R~l  2 2 

COS0  t cos  0 (cos  0 - Rsin  0)cosRA 


V ■ - sS5co8^0^  sin0  slnRX 
10 


Description  of  Variables 


AL 

ARE2 

CRL 

CT 

CTR 

CTRM 

CT2 

CTR2 

Cl 

C12 

C2 

DSJ 

DTHJ 

DXJ 

DYJ 

DZJ 

I 

J 

J1 


aj  + (l-l)AXj 
r2 

cos(RX} 

CO80J 

COS0j  (cos0j) 

. « vR-1 

(cos0j) 

COS^0j 


R-1 


cos 


2R 
’0J 

a^l/10 

w 

Proportional  to  AV  of  KH. 
A0j  - 2a 
a cos0jAXj 
a A0J 

■ “tH,  - V>, 

Zonal  Index. 

Meridional  Index. 

Loop. 
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J2  Loop. 

K Sigma  Index. 

NBJ  Temporary  variable  holds  number  of  cells  in  band  j . 

NT  Main  loop  index. 

PI  2R^  - R - 2 

P2  R^  + 2R  + 2 

P3  (R+1)^ 

RPl  R + 1 

RP1RP2  (R+l)(R+2) 

RP2  R + 2 

SRL  sin(RX) 

ST  3in0j 

THJ  Temporary  variable  to  hold  9 j . 

TR2  2r2 

UN  Returned  value  of  unit  function. 

VEL  Magnitude  of  velocity  vector. 

Helpful  Comments 

EQUIVALENCE  ((ISER,  BLK(l) , (ITS,  BLX(2)) 

ISER  - NSER  $ ITS  - NTS  $ BLK(3)  - APT (GOBALb.  13). 

Use  of  the  equivalence  avoids  BLK(l)  - NSER,  which  would  initate 
a conversion  to  floating  point.  This  3 word  identifying  label  appears  as 
the  first  record  of  the  data  file. 

ta-K2$K2«Kl->-l$  VMAX  - 0 $ KMAX  « IMAX  « JMAX  « KMIN  - JMIN  - 0 
These  values  will  emerge  as  the  indices  of  the  locations  of  the 
maximum  advecting  velocity  in  the  initial  data,  and  the  minimum  space  incre- 


ment in  the  grid. 
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J2  - 1 $ DO  10  NT  - 1,  NBLK  $ J1  - J1  + 1 $ JLI  JLI(NT)  $ LT  - 0 


:G0BAL(|>.24 


1,  NBLK 


BLK(LT) 

Reinitialize  LT  ~ 0 


ISTARl 


UN  - UNIT(ISTARl)  $ If  (UN.GE.O)  CALL  PLAY  (ISTARl.  20.  NT.  UN 


Note  that  the  second  parameter  (20),  In  the  call  to  DLAY,  locates 


the  origin  of  this  call  uniquely  within  the  outermost  DO  loop,  and  the  third 


parameter  (NT),  which  Is  the  control  Index  of  the  loop,  will  determine  on 


which  pass  through  the  loop  the  error  occurred 


Buffer  out  (ISTARl,  1)  (BLK(  ),  BLK(  ))  $ »IN  - UNIT(ISTARl)  $ LT  - 0 

VVV\A/\/VVVV\/V\AA/>rf 


The  underlined  statement  has  a history.  The  program  had  been 


this  case  there  are  (2  + 3 x 1)  x 32  x 16  ■ 2560  data  words  In  the  Initial 


data  file,  and  the  I/O  size  was  6528,  only  one  block  was  required  on  the 


ISTARl  file,  and  the  loop  Indexed  NT  was  passed  through  but  once.  However 


when  a finer  grid,  (32  x 64)  was  tried,  bizarre  results  were  obtained.  A 


10240  data  words,  and,  with 


an  I/O  size  of  6528,  this  id.ll  require  2 blocks  of  data  (NBLK  ■ 2) 


eb«  outermost  loop  Is  passed  through  twice.  Originally,  however,  the  unit 
check  was  not  performed  at  the  end  of  the  buffer  out  In  G0BAI4*64.  Thus, 
new  data  was  generated  and  placed  on  top  of  old  data  In  BLK  before  the  buffer 
was  Initialized,  le.  the  old  data  was  destroyed  before  It  was  output.  Unit 
checks  will  not  be  made  at  certain  crucial  points  In  GEX,  to  more- efficiently 
overlap  computation  with  I/O.  But,  a lesson  here  Is  that  one  must  exercise 
care  In  attempts  to  Increase  efficiency  this  way. 


SUBROUTINE  COMPflUl.  IU2) 

ENTRY  EXPAND 

EXPAND  Is  called  from  INITE  with  the  formal  parameters  (ISTARl,  IDES)  If 
the  print  option  Is  In  effect.  In  this  case.  It  Is  necessary  to  reformat 
the  ISTARl  file  (It's  In  compressed  format),  to  prepare  It  for  acceptance 
by  REPORT.  REPORT  expects  data  for  a time  step  NTS  to  be  organized  into  JBl 
records,  plus  4 dummy  records.  The  first  record  Is  a 3 word  label,  and 


subsequent  records  are  full  bands  of  data,  output  In  arrays  P(MP,  K2,  NBJ) 
(See  KMII  p.  7).  That  Is,  the  data  for  band  J Is  output  as: 

t 


EXPAND 


Since  the  CBC-6300  scores  the  S-disiensional  array  as 
sequential  planes  ?(iI,K)T,  and  stores  each  plane  as  a l-dinensional 
array  with  the  first  subscript  varying  fastest,  the  actual  storage 


where 


3-D  Variable 


2-D  Variable 


Grid  Parameters 


I • zonal  index 


J * meridional  index 


K • vertical  index 


Data  filed  coming  out  of  EXPAND,  and  active  data  files  during  GEX  runs, 
will  also  have  3 dummy  bands  at  the  end  organized  as  follows: 


J - JBl 


J=J32 1 

3 word  dummy 
record 


.J-JB4 


's  open 


Common  Blocks 


/BUF/  , 

/UNITS/ 

/ / The  array  P(4,  3,  64)  is  used  to  organize  the  data  into  the  format 

discussed  above. 


{ 


CLX  Three  word  array  used  to  hold  label  when  connress  entry  is  called. 
In  this  case  a serial  number  check  is  performed  before  the  file  is 
output  as  a display  file. 

FI  Zonal  index  (real).  (See  XMI  p.  32) 

I Zonal  index  (integer) . 


ITAl 

^l,i,j 

(KMI) 

ITA2 

'^2,i,j 

(KMI) 

IXZN 

r(XN)bj 

^INT 

Cbj-i  + 

(b,_i/bj)(Oj-, 

u J.  i-rINT 

•j)  + bj.;^  - + .o_. 

IXES 

[:(xs)bj 

-fINT 

Lbj+i  + 

(bj^^/b.)(a.- 

-fINT 

.5)  + b.^^  - + .5^ 

IXIl 

^l.i.j 

1X12 

^2,i,j 

J Meridional  index. 

J1  Control  on  loop  expanding  data  block  NT  into  ? arrays. 

J2  Control  on  loop  expanding  data  block  NT  into  P arrays. 

K Level  index. 

ST  Loop  index  for  buffer  on  dummy  bands. 

L Index  of  JLI  array  when  compressing  data. 

LT  Index  of  array  3LI<  into  which  data  is  being  compressed,  or  out  of 

which  data  is  being  expanded. 

L2  Holds  size  of  successive  data  blocks  in  compressed  file. 

M Index  for  3-dimensional  variables. 

.NBJ  bj 
NBMl  b,_, 

NBPl  bj^.3_ 

) NHOLD  Equivalent  to  CLX(2).  .Avoids  ITS  ■ CLS(2)  in  compress. 

NT  Index  for  loop  on  compressed  data  blocks. 

a / 


NUN 

Vjj  (KMI) 

NOS 

Vg  (KMI) 

SBl 

bj.i/bj 

■ 

RB2 

bj^l/bj 

XAl 

bj_i/bj  (aj.  5)  + + .5 

TA2 

^j+l^^j  '*■  '’j+1  ” ®j+l 

UN 

Returned  value  of  unit  function. 

NN 

’'axjSi 

WS 

XN 

i(bj_i/bj)  + (bj_i/bj)(Oj_  3)  + bj_^ 

- “j-i  + -5 

ZS 

l(bjVb^)  + (bj^l/b^)(aj_5)-.b^^l 

- V" 

Helpful  Comments 

ENTRY 

EXPAND  (GCE. 32) 

Calls  to  COMP  Xflll  be  aaalyzed  in  the  natural  sequence  they  will  occur. 

C EXPAND  FILE  ON  lUl  BUFFER  TO  IU2 (GCE.  33) 

Note  that  the  £lle  Is  assumed  positioned  by  the  calling  program. 

NTS  ■ ITS  $ APT  - BLK(3)  $ NSER  - ISER (GCE. 37) 

Save  3 word  label  off  this  data  file.  In  a GEX  run,  with  restart  data 
in  two  time  steps  ISTARl  and  ISTAR2,  the  first  file  expanded  is  step  N-1, 
the  next  Is  step  N.  In  that  case  the  statement  above  leaves  NTS  ■ N on 
return  to  GEX  (as  It  should!),  and  the  next  time  step  prognosticated  will 
be  Nfl. 

J2  - 1 $ DO  200  NT  - 1,  NBLK  $ LT  - 0 
LZ  - LSZOn)  

The  descriptive  variables  for  the  compressed  file  are  here  used  to 

organize  chat  file  into  an  expanded  file,  of  1 label  record,  followed  by 
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JB  data  records , followed  by  4 dummy  records.  The  LSZ  array  carries  the  size 
of  each  of  the  NBLK  compressed  data  blocks. 

The  loop  begins  with  the  compressed  data  file  positioned  after  the  label, 
and  before  the  first  data  block. 


The  logic  of  the  loops  Is: 


K2  are  fixed  in  P (MX , K2 , NBJ)  ,•  but,  NBJ  nay  vary  band  by  band.  MX  is  used 
rather  than  MP,  since  the  prognostic  variables  occupy  fixed  positions  within 
each  planar  segment  of  P. 


TAl  - RBI  * (DEV(J)-.5)  + FNB(J-l)  - DEV(J-l)  + .5 

TA2  - RB2  * (DEV(J)-.5)  FMB(J-H)  - DEV(J-H)  + .5 (GCE.51) 

That  is  TAl  - (Oj.  5)  bj_j^/bj  + bj-l  “ Oj-l  + *5 

TA2  ■ * ^J+1  ” ®j+l 

IXEN  " m(J)  * RBI  + TAl  $ IXES  “ FNB(J)  * RB2  + TA2 (GCE.52) 

That  is  IXEN  ■ + (Oj-.s)  ■*"  ^j-1  “ ®j”l 

IXES  ■ bj^j^  + (aj_  j)  + .5 
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Notice  the  conversion  to  Integer. 


Now  - l(bj_i/bj)  + *^0 


(See  KMl  3.5.3). 


Therefore  IXEN  - 


■ l*N-l  \ 


That  Is,  these  values  are  Initialized  to  that  corresponding  to  the  last 
cell  In  the  band  J with  Index  bj , before  entering  the  1 loop. 

ITA2  ■ MOD (IXEN, NBMl)  + 1 $ 1X12  - MOD(IXES.NBPl)  + 1 

. rr^-TiimMOD('»j-i)  . , 


+ 1.  (KMI  3.5.4) 


FI  - I $ XN  - FI  * RBI  + TAl  $ XS  - PI  * RB2  + TA2 


These  are  the  forms  for  Xj^,  Xg  given  above  for  cell  1. 
P(3.K1.I)  - WN  * AM0D(XN,1.)  $ P(3,K2.I)  - WS  * AM0D(XS,1. 


I r M0D(1) 

That  is  P(3,K1,I)  - Xg  - M 


P(3,1C2.I)  - X|  - W^xf+i 
(See  KMI  3.72) 

ITAl  ■ ITA2  $ IXIl  - 1X12  (GCE.57) 

That  Is 

(See  KMI  3.3.11). 

On  the  first  pass  through  the  loop,  for  example,  thin  statement  sets 
the  left  northern  boundary  Index  of  the  first  cell  on  the  band  equal 

to  the  right  northern  boundary  Index  of  the  last  cell  on  the  band 
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j 

■ I 

SDBROUTINE  NU(I1.I2.NB.NUT) 

* This  routine  returns  the  coded  value  of  v^,  Vg  for  the  cell  1 on  the 

band  J.  Notice  that  since  the  grid  If  organized  In  non-overlapping  vertical 

colums,  these  values  will  serve  for  all  k levels. 

Description  of  Variables 

n Westernmost  contiguous  cell 
For  north  boundary 

For  south  boundary  (5£  j j+1^ 

12  Easternmost  contiguous  cell 

For  north  boundary  (^2,i,J,j-l) 

For  south  boundary  ^C2,i,j,j+l) 


NB 

Number  of  cells  on  contiguous  band 

For  north  boundary  bj.^ 

For  south  boundary  bj.|.]^ 

NUM 

Temporary  variable  to  hold 

coded  value  of 

NUT 

Coded  value  of  v.  See  description  ahead. 

NUl 

’’2  " ^C2  " 

■ ' 

NU2 

bj^l  + n2  - (ni+1)  or  bj^j^ 

+ C2  - (Ci+1) 

NUT 

Value  NUT 

V ■ 0 Easternmost  and  westernmost  cells 
are  contiguous 

V -1  Eastermost  and  westernmost  cells 
coincide 

V - Hj  - ni  or  52  - Cl 

bj.i  + n2  - (ni+1)  > 0 (bj+2  + C2  - (^1+1)  ^ 0) 

bj.i  + n2  - (n^+l)  > 0 (bj^.2  + 5i  - (Ci+1)  1 0) 

”1  **  **J-1  ^^1  ^j+1^  _ 

n2  S 1 CC2  i 1) 
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6 + '^2  “ ^ ° ^^J+1  ^2  ■ ^ 

n^  < bj-i  (Cl  < bj+i) 

112  > 1 ^^2  ^ 


Selected  data  from  file  lU,  which  has  just  been  expanded,  is  reported 
on  the  community  output  file.  REPORT  is  called  from  DISPLAY  with  MR  - 1 or 
MR  • -1,  and  from  INITE  with  MR  ■ 2. 

Common  Blocks 
/BUF/ 


/UNITS/ 

/AC/ 

/ / 
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Variable  Deacrlptlon 


CT  cos(0j)AXjAej 

I Zonal  Index, 

in  Expanded  file. 

J Meridional  Ixdex. 

JB3  JBl  + 1 

JH  J - 1.  Current  value  of  J being  processed. 

JF  Data  index  of  band  1 JF  JB. 

JTS  Equlvalenced  to  PI (2).  Used  to  hold  NTS  from  label. 

K Vertical  Index. 

M Variable  Index. 

M - 1 u 

M - 2 V 

M - 3 T 

M - 4 r 

MM  Index  for  output  array. 

MPl  MP  + 1 

MF2  MP  + 2 

MR  -1  Just  report  TE  and  TM 

0 Just  report  label  and  real  time  clock 

1 Called  from  DISPLY 

2 Called  from  INITE 

NBJ  . bj 

PE  Partial  energy.  Energy  in  band  J. 

PM  Partial  mass.  Mass  In  band  j. 

TE  Total  energy. 

TM  Total  mass. 
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Th«  logic  of  the  loops  Is: 


K - 1,  KV 

PE  - 0,  PM  - 0 

1 - 1,  NBJ 

I M - 1,  MP 

P3(M)  - Pl(M,K,i) 

I 12 

P3(MP1)  - 
P3(MP2)  - 
PE  ■ PE  + 

PM  - PM  + 

PE  - PE  + 

IF  (HR  ± 0)  GO  TO  11 
WRITE  P3 
^11 


Pld.Kl.D- 
P3(MP1)  - 8000 
P*(u2  + 

P* 

go  (P*)' 


v2) 


PM  ■ PM*cos0jAXjA9 


J 


PE  ~ PE*cosd jAAj Ad4 
WRITE  PE,  PM 
-10 


3-D  prognostic  variables 
for  cell  (i,J,_k) 

• P^  for  (i.j.k) 


Helpful  Conments 

EQUIVALENCE  (Pl(2).  JTS) (GRPORT.20) 

JTS  holds  time  step  number  after  label  is  buffered  into  PI. 

REWIND  lU 

After  leaving  expand  the  file  lU  is  positioned  at  the  end  of  the  last 
dummy  band  JBS. 

WRITE  (6.  200)  ITS,  BLK(3).  RTJ  fGRPORT.29) 

Before  leaving  EXPAND,  the  label  was  returned  to  BLK(1)~BLK(3) . 

39  BUFFER  IN  (lU.l)  (PKl.l.U  .P1(MX.K2.NB5)) 

Only  the  3 word  label  is  buffered. 

J - 2 

NBJ  - NB(5) (GRPORT.49) 

J - 2 starts  the  first  data  band,  with  NBJ  cells. 

BUFFER  IN  (IU,1)  (P2 (1.1.1).  P2CMX.K2.NBJl) (GRP0RT.56) 

Notice  that  there  is  not  a unit  check,  immediately  after  this  buffer. 

In  this  way  the  data  from  the  band  (PI),  can  be  manipulated  simultaneously 
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vhlle  the  buffer  for  the  J+1®*^  band  (P2)  is  being  initialized. 

JM  - J - 1 

HBJM  - NB(JM) (GRP0RT.59) 

JM  is  now  the  meridional  index  of  the  band  in  PI. 

PE  - 0.  $ PM  ■ 0. (GBPORT.ei) 

Energies  and  masses  are  reported  for  each  band  at  each  fixed  level. 
Here  partial  energy  and  mass  are  reinitialized  before  reporting  the  next 
level. 

P3(M)  - Pim.K.I)  (GRPORT.64) 

P3  la  used  only  to  hold  the  selected  variables  to  be  printed. 

P3aiP2)  ■ P3fMPl)  - 8000. 

A convenience  for  scanning  the  output. 

PE  - PE  P3(MP1)  * (P3(l)**2  P3(2)**2)  $ PM  - PM  + P3(MP1) 

PE  ■ PE  GO  * P3(MP1)**2  (GRPORT.69) 

Energy  « Kinetic  Energy  Potential  Energy 

- E P*  { (u2  + v2)  + g^jP*} 
ij.k 

In  the  baratroplc  case  P*  ■ H ■ height  of  free  surface 

Mass  ■ E P.  where  (P*)<  4 “ H.  . - AMass  above  (l,j)  (See  4.6  KH) 
IJ.k* 

JP  - JM  - 1 (GRPOriT.  70) 

JP  is  the  index  of  the  present  band  in  the  data  bands  (1  ^ JP  JB 
whereas  1 £ J JB5) . 

PM«PM*CT$PE«PE*CT$TE-TE-»-PE$TM-TM-hPM 

■ <“>j,k  “l.! 

where  (4M)  j is  mass  per  unit  area 

(dE)4  is  energy  per  unit  area. 


60 


■ cos0jAXjA0j  weights  the  cell  areas  for  band  J. 


WRITE  (IPRT.  700)  NSER,  NTS,  JP,  PE,  PM 


NTS  was  saved  as  the  value  of  JTS.  Partial  energies  and  masses  are 


reported  after  the  prognostic  variables  for  each  band,  and  for  each  level. 


IF  (J  - JBl)  15.  16,  16 


(GRP0RT.81) 


When  J » JB3  * JBl  + 1 =•  JB2  all  data  bands  have  been  reported. 


17  BUFFER  IN  (lU.l)  (Pl(l.l.l).  PI CMX.K2 .NBJ) ) 


GRPORT.85) 


Notice  again  no  unit  check  immediately  after  the  buffer.  The  data  from 


the  J+1®  band  (P2)  can  be  manipulated  (simultaneously)  while  the  buffer  for 


the  J+2^  band  (PI)  is  being  initialized. 


JM  - J - 1 
NBJM  - NB(J-l) 


(GRP0RT.87) 


JM  is  the  meridional  index  of  the  band  in  ?2. 


APPENDIX  B 


"Global  Experiment" 


GEX  Integrates  the  variables  in  time  and  saves  user  specified  time 


steps  for  later  display. 


Characteristics 


Core:  70000B 


Time : Open 


Options:  Initial  data  or  restart  data 

Choice  of  time  integration  schemes 
Choice  of  time  steps  to  be  saved 


Permanent  Files:  Display  file 
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Organization  of  Core  During  a Rim 


Named  Common 

/BUT/ 

101 

14600 

/UNITS/ 

14701 

15 

/AC/ 

14716 

2260 

Programs 

GEX 

17176 

4766 

REPORT 

mi 

/FRC/ 

/LTC/ 

/CEH/ 

STIS(^ 

/FG/ 

/BC/ 

ICEF/ 

CBUFN  35632  212 


SYSTEM 


r 


Blank  Common  / / 54746  12000 

Common  Blocks 

/AC/  As  In  miTE. 

/BUP/  Used  as  buffer  area  for  output  (Input)  of  display  files. 

/UNITS/  As  la  INITE. 

/ / Each  PJ  array  holds  1 latitude  band  In  expanded  format.  There 

are  a ^ul■x^^ml1n  of  12  bands  needed  In  core  (for  TIS  Euler  trapezoidal) . 

Files 

INPUT (5) — Community  Input  File 

Output (6) — community  Output  Pile 

INIT(7) — Private  Binary  File  of  Grid  Information 

ISTAR1(8) — Initial  data  file  created  In  OBAL  (ts  ■ 0)  or  restart  data  file 
(ts  - N - 1) 

ISTAR2(0,  8,  or  other)— -Not  used  (0)  or  restart  data  file  sequentially 

following  on  ISTARl  file  (8) , or  found  on  another 
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INPUT  OUTPUT  CHART 


IDES (9) — Display  file  of  compressed  data  stored  as  a permanent  file 
INI  (10,  11,  12) — Active  data  file  holding  ts  .V. 

IN2(10,  11,  12) — Active  data  file  holding  ts  N-1. 

IOT(10,  11,  12) — Active  data  file  holding  ts  N+1. 

Description  of  Variables 

FTI  Fixed  time  of  integration.  Read  in  on  PASAM  card,  it  determines  the 
value  of  AtCDT)  for  this  run. 

ID  70  character  identification  for  this  run. 

IDDY  Holds  value  of  IMl  during  process  of  rotating  file  names. 

IDY  Holds  value  of  INI  during  process  of  interchanging  file  names  INI  and 

IN2  holding  restart  data. 

INPUT  Community  input  file. 

INI  Active  data  file. 

IN2  .-Active  data  file. 

IRE  Nor  presently  used. 

L Index  for  output  of  identifier  array  ID. 

U Returned  value  of  unit  function. 

UN  Returned  value  of  unit  function. 


! 

Helpful  Counments 

, 29  CALL  EXPAND  (ISTARl.  INI) (GGEX.  33) 

Tine  step  N - 1 is  organized  In  expanded  format  of  JBS  records  (label, 
JB.data  bands,  and  4 dummies).  Label  INI  Is  temporarily  attached  to  N - 1 
file. 

IF  (ISTAR2  .EQ.  0)  GO  TO  30  $ IF  (ISTARl  .NE.  ISTAR2)  GO  TO  26 

The  default  value  of  ISTAR2  Is  zero  In  the  data  statement.  If  this 

Is  not  an  Initial  data  run,  ISTAR2  should  be  set  to  8 If  time  step  N Is  on 

the  sane  file  as  time  step  N - 1,  and  to  LFN,  If  time  step  N resides  on  a 

distinct  file  LFN.  If  this  Is  an  Initial  data  run  an  Euler  step  Is  done  first. 

BUFFER  IN  (ISTARl,  1)  (BLK(l),  BLK(3))  $ UN  - UNIT(ISTARl) 

IF  (U  .NE.  0)  GO  TO  750  (GGEX.  36) . 

Skip  over  EOF  If  both  N and  N ~ 1 on  same  file.  If  unit  check  does 

not  return  zero,  write  error  condition  and  stop. 

26  CALL  EXPAND  (ISTAR2.  IN2)  $ IDY  - INI  $ INI  - IN2  $ IN2  - IDY  $ GO  TO  35 

For  restart  data  option  time  step  N Is  organized  In  expanded  format  of 

JBS  records.  Label  IN2  Is  temporarily  attached  to  N file.  Since  N - 1 step 
was  created  first  In  a preceedlng  GEX  run.  It  will  be  the  first  off  the  dis- 
play file.  At  the  point  above,  the  file  names  are  switched  so  that  IN2  holds 
ts.  N - 1,  and  INI  holds  ts.  N,  as  we  wish.  The  branch  to  35  skips  the  Euler 
step  and  goes  directly  to  a TIS  step. 

30  CALL  STIS<^ 

CALL  ROTATE  (TIS(») 

For  Initial  data  option  these  statements  Initialize  an  Euler  step  and 

put  the  generated  data  In  expanded  format  on  file  lOT.  Note  that  TIS^  and 

» 

I TISl  are  subroutine  names  In  an  external  statement. 


L 


I 


34  IDDY 

- INI  $ 

INI  - lOT 

$ lOT  - 

IN2  $ IN2  - IDDY 

To 

avoid  the  actual 

transfer 

of  data,  the  file  names  are  rotated. 

this  way 

(12) 

IN2 

N - 1 

(11) 

INI 

N 

* 

(10) 

lOT 

N + 1 

becomes 

(11) 

IN2 

N 

• 

(10) 

INI 

N + 1 

(12) 

lOT 

N - 1 

and  the  stage  Is  set  for  the  creation  of  ts.  N + 2 in  lOT. 


35  CALL  STISl 

CALL  ROTATE (TISl) 


A TIS  step  Is  initialized  and  the  generated  data,  in  expanded  format, 
is  put  on  file  lOT.  The  nature  of  the  TIS  (Leapfrog,  Adams-Bashforth,  Leap- 
frog Trapezoidal,  etc.)  Is  determined  by  the  update  * COMPILE  card.  All 
modules  share  the  same  subroutine  names,  but  have  unique  deck  names,  enabling 
the  user  to  compile  a specific  TIS  from  his  Global  library,  and  then  use 
the  above  Fortran  calls  to  execute  that  TIS. 


IF  (NTS  .LE.  MTS)  GO  TO  34 

Note  that  the  return  to  3A  rotates  the  file  names,  before  pressing 
ahead  with  the  creation  of  the  next  ts.  In  lOT. 

SUBROUTINE  STIS(» 

Start  time  Integration  scheme  Euler.  This  routine  sets  up  data  arrays 
In  central  memory  In  preparation  for  going  Into  TIS  at  J * 2 (first  data 
band) . 
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Coaccn  Blocks 


As  previously  described 


/ / 

/BOP/ 

/DHITS/ 

/AC/ 

/LTC/  Codes  and  ocher  InforasClon,  used  In  creating  display  files  while 


Clae  Integration  scheme  Is  proceeding. 

/CEP/  Grid  Information  for  a particular  cell  (l»j(k).  This  Information  Is 
passed  to  Che  routines  which  calculate  Che  tendencies  for  that  cell  (which 
are  also  held  In  this  common  block). 

/CEB/  Holds  bj  for  a fixed  band  being  processed. 

/CEF6/  Holds  the  zonal  Index  for  a specific  column  of  cells  being  processed. 
Description  of  Variables 

COP  Coriolis  Parameter  f ■ 2Qsln6j  cos(A6j/2)  for  current  band  being  pro- 
cessed. 

COT  Metric  term  In  equation  of  motion  m ■ tanOj/a  for  current  band  being 
processed. 

D3  Integral  of  3-dlmenslonal  divergence  for  the  A prognostic  variables  at 
25  possible  levels. 

1 Holds  zonal  Index  for  current  cell  column  being  processed  by  RHS. 


ICE  Code  to  locate  origin  of  call  to  buffer  error  subroutine,  within  this 
subroutine.  These  codes  are  unique  to  each  routine  (See  appendix  for 
BERROR  codes). 

IDE  Code  to  Indicate  If  the  time  step  currently  being  processed  Is  a display 
fUe. 

IDX  * 0 Not  a display  file 
IDX  - 1 Display  file 

IM  Zonal  Index  of  cell  Immediately  west  of  current  cell. 


lOPS  Code  for  sense  switch  I 1. 

lOPS  " 1 Operator  stop  Indicated 

ICPS  ■ 2 Processing  continues. 


IP  Zonal  index  of  cell  immediately  east  of  current  cell. 

ITA  Holds  n][i  n2  for  current  cell  column — ITA(l)  ■ ITA(2) 


IHN 

ZXZ 

JT 

LB 

LT 

NBJ 

NBM 

NBP 

NBl 

NB2 

NB3 

NUN 

NUS 

RAP 

RT 

RTE 

WEVJ 

WNJ 


Identifies  file  on  which  buffer  error  occurred — passed  to  BERROR 
subroutine. 


Holds 


^1*  ^2  current  cell  column — IZI(l)  “ 1X1(2)  • 


Identifies  call  to  BERROR  routine  as  having  originated  in  STIS  sub- 
routine (JT  ■ 0). 

Index  for  JLI  array.  If  current  ts.  is  a display  file  (IDX  - 1),  then 
LB  is  incremented  each  time  a block  of  data  goes  out  as  a record  on 
the  display  file. 

Index  for  BLK  array.  If  current  ts.  is  a display  file  (IDX  ■ 1),  then 
LT  is  Incremented  as  BLK  is  being  filled  with  compressed  data  to  be 
output  to  display  file. 

Number  of  cells  in  current  band  bj . 

Number  of  cells  on  band  north  of  current  band 

Number  of  cells  on  band  south  of  current  band 

Holds  bj  for  first  dummy  band,  NBl  ■ 1. 

Holds  bj  for  first  data  band,  J - 2. 

Holds  bj  for  second  data  band,  J > 3. 

Holds  coded  for  current  cell  column. 

Holds  coded  Vg  for  current  cell  column. 

3P. 


P*  tendency. 


3t 


for  current  cell  column. 


Real  running  time.  If  RT  ^ RTL  - 10,  end  processing  is  initiated. 
Tendencies  of  3-D  prognostic  variables  for  each  of  the  KV  levels  indexed 

(I.J). 

East-west  weight  for  current  band  (KMI). 

North  weight  ^ for  current  band  (KMI) , where  Al  is  the 

2 

length  of  the  north  coincident  boundary  segment,  and,  AAj  ■ 2a  AXj 
cos0jSln(^®j/2)  is  constant  for  each  band. 
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f WHSJ  - Wg,  where  ■ ZW^,  and  Wg  ■ EW^.  Value  appears  in  Gg  operator  (KH). 

[ US 

' > HSJ  South  weight  ~ dl/AAj,  for  current  band, 

i HXN  Holds  for  current  cell  column,  WXN(l)  ■ UXN(2)  ■ 

[ HXS  Holds  W^2  current  cell  column,  vncs(l)  - WXS(2)  - 

Helpful  Comments  ‘ 

LT  ■ 0 $ IF  (ISS  .EQ.  0)  GO  TO  1011  $ WRITE  (IPRT.  1012)  NTS  $ GO  TO  1007 

On  first  entering  STIS^,  ISS  has  been  preset  zero  In  GEX.  On  subsequent 
^ entries  ISS  may  have  been  set  to  1,  If  TIS^  has  detected  a negative  value  In 

the  scalar  field. 

1007  D - UNITflNl)  ■»•  UNIT(IN2)  $ REWIND  INI  $ REWIND  IN2  $ CALL  CC«g(IM2.1DES) 
If  any  one  of  the  stop  conditions  exist  then  INI  and  IN2  are  compressed 
on  the  display  file.  Since  the  file  names  were  rotated  in  GEX  before  the 
entry  to  STIS^,  INI  contains  the  last  time  step  created,  and  IN2  the  previous 

one.  Note  that  the  display  file  counter  Is  incremented  two.  The  subroutine 

COMP  also  has  the  entry  expand. 

1001  NTS  - NTS  + 1 $ ADT  - ADT  + DT  $ IF  (NTS  .NE.  MTS)  GO  TO  1006  $ 

IF  (IPX  .NE.  0 GO  TO  1009 

None  of  the  conditions  for  an  immediate  stop  has  been  met.  Since  the 

next  time  step  will  be  created,  NTS  Is  Incremented.  A test  Is  made,  to  see 

If  this  will  be  the  last  time  step.  If  so,  both  the  ts.  NTS  and  NTS  - 1 must 

I be  saved  on  display,  though  neither  may  be  a regular  display  file.  This 

, la  necessary,  since  we  need  two  consecutive  steps  for  restart.  Thus,  a test 

Is  made  to  see  If  NTS  was  already  saved  (IDX  "I).  If  not,  INI  Is  compressed 

j on  display  and  lOT  Is  coded  to  be  displayed  (IDX  ■ 1).  Otherwise,  lOT  Is 

I coded  to  be  displayed,  and  Its  label  sent  out  on  IDES  file. 

I i 
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j 

I 

1009  ISER  » NSER  $ ITS  « NTS  $ BLKQ)  - APT  $ IPX  - 1 $ CALL  PBUFT(IDES.95.0.a.BLK3^ 
Codes  ts.  NTS  to  be  displayed,  and  puts  the  label  on  the  file  IDES.  The 
code  95,0  Is  used  In  error  processing  by  the  BERROR  routine.  Note  again, 

ISER  - NSER  avoids  BLK(l)  - NSER. 

IF  (UNIT(INl))  1.2.3 

Note,  here  and  below,  that  each  type  of  error  and  location  where  It 
occurred  has  a unique  ICE  code. 

BUFFER  IN  (INl.l)  (PSa.l.l).  P3CNX.K2.NB2))  (GSTISE.9) 

The  label  off  file  INI  Is  placed  In  P3  array. 

CALL  PH1(P3.  P3.  PHI) (GSTISE.  19) 

This  call  to  the  geopotentlal  routine  Is  made  with  the  dummy  label  In  P3. 

PI (3)  “ BLK(3) (GSTISE.  25) 

The  label  for  the  cxirrent  time  step  Is  placed  In  array  PI.  Then  on 
leaving  STIS(^  we  have: 


PI 

P2 

P3 

P4 

P5 

Label  for 
current  ts 

NTS 

Arbitrary 

Label  for  ts 

NTS  - 1 

data  band 
for  ts  NTS  - 1 

J - 2 

2nd 

for 

J - 

data  band 
ts  NTS  - 1 
3 

PHI  PH2 

Arbitrary  4 field 

for  1^^  data  band 
at  all  K7  levels 


SUBROUTINE  COMPCIUl.  IU2) 

Compress  file  lUl  Into  file  IU2.  See  entry  EXPAND  for  a description  of 
the  data  format  In  the  expanded  file  lUl;  see  INITE  for  a description  of  the 
data  format  In  the  compressed  file  IU2. 

i 

* 

? 
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/UNITS/ 


/BUT/ 

/ / 


Variable  Description 


Helpful  Comments 


ISER  - NSER  $ IF  (BLK(l)  .NE.  CLKa)LGO  TO  201 


A serial  number  check  on  the  data. 


(GCE.12^ 


BUFFER  OUT  (IU2.1)  BLK(l).  BLK(3 


GCE.15 


Label  goes  out  on  IU2,  exactly  as  it  came  In  on  INI. 

The  loop  begins  with  the  expanded  file  lUl  positioned  after  the  label, 
and  before  the  first  data  band.  The  logic  of  the  loop  Is: 


LT  Is  reset  to  zero  whenever  BLK  is  to  be  refilled  with  compressed  data. 

L keeps  count  of  blocks  being  created,  the  total  number  of  which  was  predeter- 
mined in  the  IMIT  run. 


IF  (J  .LT.  JLKL))  GO  TO  100 


(GCE.26 


Go  to  100  if  this  band  is  not  the  last  latitude  in  a compressed  block; 
otherwise,  output  the  data,  since  the  next  band  in  will  exceed  the  I/O  size. 


SUBROUTINE  CBUFN  (lU,  I.  J.  U,  A.  N 


ENTRY  PBUFN 
ENTRY  CBUFT 
ENTRY  PBUFT 

N words  of  data  are  buffered  in  (subroutine  name  ends  in  N) , from  lU; 


or,  buffered  out  (subroutine  name  ends  in  T)  to  lU.  If  prefix  is  C,  unit 


checks  are  done  on  entry  and  on  exit,  if  prefix  is  P,  unit  checks  are  done 
only  on  entry.  CODE  is  passed  to  BERROR  if  buffer  error  occurs. 


SUBROUTINE  PHICPJ,  PJS,  PHJ] 


Returns  geopotential  scalar  field  for  band  j at  all  KV  levels.  The 
value  of  geopotential  is  given  at  the  top  and  bottom  of  a box.  Subroutine 
is  now  programmed  for  1 level.  In  the  general  case  PHJ  should  be  dimensioned 
PHJ(m,NBJX). 

Common  Blocks 
/BUF/ 

/UNITS/ 


/CEH/ 


Helpful  Comments 


PHJ(KV1,I)  - GO  * PJ(1,K1,I) 

1 PHJ(l.I)  « PHJ(KVl.I)  (GPHIBR. 10] 


„ _ AifO 

From  ^ ~2  * 


*5(®lc+i5  + Ok-ls^ 


we  have 


♦k-ls  ■ ♦kWs  * 


‘•‘“kW,  * Vli> 


Since  we  carry  field  ■ '^kv+Js^  ’ determine  the  (p  field 

Iteratively  through  the  top  of  the  atmosphere  (t^^) . In  the  special  case 
here  programmed,  Tj^  = 0 in  space  and  time,  so  that 

♦k-»s  * 

• Vh  — PHJ(KV1,I)  - - goH 

n la  Che  height  of  Che  free  surface,  and  is  carried  in  PJ(1,R1,I). 

ei  ee . PU  ' 1 . 1 : • >,  • goH 
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Note:  calls  to  FHI(PJ,FJS,PHJ)  with  dunmy  bands  as  PJ,  return  data  of 
no  significance  In  PHJ. 

SUBROUTINE  ROTATE  (TIS) 

Rotation  Subprogram.  This  routine  rotates  the  data  arrays  from  the 
north  pole  to  the  south  pole  In  the  process  of  updating  the  prognostic 
variables  to  the  next  time  step,  with  a chosen  procedure  subroutine  TIS. 

The  logical  rotation  of  the  data  avoids  the  actual  transfer  of  data  In  main 
memory.  When  the  data  units  are  left  rotated,  they  assume  the  proper  con- 
figuration for  updating  the  next  latitude  band  j + 1.  A maximum  of  12 
prognostic  arrays  and  8 diagnostic  arrays  Is  necessary  to  accomodate  the 
most  complex  of  the  IIS’s. 

Common  Blocks 

/ / Holds  arrays  which  are  being  rotated. 

/BUF/ 

/UNITS/ 

/AC/ 

/BC/  Holds  Index  of  current  latitude  band. 

Variable  Description 
No  new  variables  Introduced. 

Helpful  Comments 

J - 1 

36  J - J 1 (GROTATE.28) 

The  first  band  processed  will  be  Index  J * 2 (the  first  data  band).  A 
return  to  36  will  be  made  If  a rotation  back  to  the  original  order  Is  called 
for,  and  not  all  the  data  bands  have  been  processed.  Since  the  least  common 
multiple  Q2,  8]*  24,  this  will  occur  after  23  rotations  (at  statement 


GROTATE. 147). 


On  the  first  pass,  performs  time  Integration  on  the  first  data  band 
(J  * 2).  The  P and  PH  arrays  contain  that  Information  Indicated  at  the 


end  of  STIS. 

IP  (J  - JBl)  1.  2.  2 

Tests  If  all  data  bands  have  been  processed.  If  not,  latitude  Index 
J Is  Incremented,  another  rotation  of  data  arrays  Is  performed,  and  the 
next  band  Is  processed.  If  all  bands  have  been  processed,  the  subroutine 
ENTIS  Is  called  to  end  processing  on  this  time  step,  by  putting  out  the 
last  updated  data  band,  which  should  reside  In  array  P2,  and  also  putting 
out  the  dummy  bands.  Then  a return  to  GEX  Is  called. 

SUBROUTINE  TIS(t 

Euler  time  Integration  scheme.  This  routine  updates  the  Information 
of  latitude  band  J to  time  N. 

Common  Blocks 
/LTC/ 

/ / 

/BUF/ 

/UNITS/ 

/AC/ 

/FG/  Holds  parameters  for  line  Integration  around  cell,  for  all  KV  cells 
of  current  (l,j). 

/BC/ 

/CEF/ 

/CEFG/ 
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Variable  Description 


DAI  ■ P4(1,K2,I)  Contains  both  n2  and  code  for  current  index 

j.  and  - DAI  - [dA^^”^ 

DA2  ■ P4(2,K2,I)  Contains  both  ^2  ^Jnd  Vg  code  for  current  index  (i,j). 

[PAU^**^  - C2,i,j  Vg  - DA2  - [pA2]^^ 
lUN  Holds  file  name,  on  which  buffer  error  occurred. 

R Level  index. 

L Zonal  index  of  cells  contiguous  to  current  cell  (i,J)  as  line  inte- 

gration proceeds  along  cell  boundary  in  RHS. 

M Index  of  3-D  prognostic  variable  in  P array. 

MB  Code  passed  by  RHS  to  line  integration  routines. 

MB  * 1 Integrate  on  east-west  boundary 

MB  • 2 Integrate  on  north-south  boundary 

M Do  Ickop  index  for  putting  out  display  file. 

NBdMl  ^J~l* 

NBJrZ 

PH  Array  holds  ) for  all  KVl  level  boundaries, 

2 k-%  k-*S 

and  current  index  (i,j).  Returned  from  subroutine  PSP,  which  is 
called  by  RHS. 

PST  N - 1 value  of  P,^  for  current  (i,j). 

PSTl  Updated  (N)  value  of  P*  for  current  index  (i,j). 

FT  Array  holds  — for  all  KV  levels  and  current  index 

2 2 

(1,J).  PT  is  returned  from  subroutine  PST,  which  is  called  by  RHS. 

PY  Array  holds  ^ ^ *^1  for  all  KV  levels  and  current  index 

2 2 

(1,J).  Here  ij»  is  one  of  the  3-D  prognostic  variables.  PY  is  re- 
turned from  subroutine  FSY,  which  is  called  by  RHS. 

YP  Array  holds  ^ ^.1  for  all  KV  levels,  for  current 

index  (l,j),  and  each  possible  3-D  variable  i];.  Y is  u along  east- 
west  boundary,  and  v along  the  north-south  boundary.  YF  is  returned 
from  subroutine  YSL,  which  is  called  by  RHS. 
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P*o  + *0 

Array  holds  — for  all  KV  levels  for  current  index 

2 2 

(l,j).  Y is  u along  east-west  boundary,  and  v along  north-south 
boundary.  Y1  Is  returned  from  subroutine  YSL,  which  Is  called  by 
RHS. 


ENTER 


»6(MX,  K2, 
NBJP2) 


P1(MX,  K2, 

J 

NBJMl) 

Compute! 


— I-1,NJB 


Set  grid 
for  Values 
for  band  J 


Call  RHS  to  compute 
tendencies  at  KV 
levels  of  Index  (i,j) 


RETURN 


this  a dlspla 


Place  N 
data  on  10 


Place  updated 
variables  In 
array  P2 


Helpful  Connnents 

NBJP2  - NB(J+2) 
NBJMl  - NB(J-I) 


(GTISE.66) 


On  first  call  from  ROTATE,  J - 2 (first  data  band) . The  J - 1 is  the 
index  of  the  dummy  band,  and  J + 1 is  the  index  of  the  anticipatory  band.  On 
subsequent  calls,  J is  the  current  band  being  processed,  J - 1 and  J + 1 are 
needed  for  an  Euler  step,  and  J + 2 is  brought  in  as  anticipatory  data. 

I BOEFER  INCINl,  1) (P6(l.l.l) .P6(MX.K2.NBJP2))  (GTISE.69) 

No  unit  check  after  the  buffer,  to  allow  overlap  of  I/O  on  this  band 
(whose  data  will  not  be  needed  until  J is  incremented  in  ROTATE) , with 
processing  in  RHS. 

4 BUFFER  OUT  (lOT.l)  (Pl(l.l.l).  P1(MX.K2.NBJM1)) 

On  first  call  from  ROTATE,  the  3 word  label,  put  in  FI  before 
leaving  STIS,  is  put  as  the  first  record  on  file  lOT.  On  subsequent  calls 
from  ROTATE,  PI  holds  the  updated  data  from  current  band  J - 1.  This  data 
was  placed  in  array  P2  during  the  last  pass  through  TIS. 

CALL  PHI(P5.P5,  PH3) 


\ Geopotential  for  band  J + 1.  On  completion  of  this  call  we  have:! 


P2 


P3 


P4 


P5 


P6 


(J-1,N) 

(J-2,N- 

1)  (J-1,N-1) 

(J,N-1) 

(J+1,N-1) 

(J+2,N-1) 

PHI 

Ph2 

PH3 

(J-1,N-1) 

(J,N-1) 

(J+1,N-1) 

The 

present 

subroutine  will  replace  the  data 

in  P2  (which  is 

no 

longer  needed) , with  the  updated  data 

for  the  current 

index  J . Then 

before 

leaving  TIS 

we  have 

the  state: 

PI 

P2 

P3 

P4 

P5 

P6 

(J-1,N) 

(J.N) 

(J-1,N-1) 

(J.N-1) 

(J+1,N-1) 

(J+2,N-1) 

PHI 

PH2 

PH3 

(J-1,N-1)  (J,N-1) 


(J+1,N-1) 
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ROTATE  then  increments  J,  and  rotates  the  array  names.  Then  on  return 


(J,N-1)  (J+1,N-1)  (J+2,N-1) 

This  is  exactly  the  beginning  state  for  the  current  index  J + 1. 
WNSJ  ■ DWNS(J) (GTISE.  87) 


COP,  COT,  WEWJ,  WSJ,  WMJ,  WNSJ  are  grid  parameters  for  the  current 
index  J,  and  are  independent  of  i,  R.  They  are  set  outside  the  loop. 
ITA(2)  - P4(1,K2,IM) 

1X1(2)  - P4(2,K2,IM) 

WXN(2)  - P3(3,K1,IM) 


WXS(2)  - P4(3.K2.IM)  (GTISE. 95) 

Initialize  n-  ,5 


Wax? 


and  W. 


IF  (I  - NBJ)  7,  8.  9 


(1-12  bj) 


7 Current  cell  has  cells  to  east  yet  to  be  processed 

8 Current  cell  is  last  cell  on  band  to  be  processed.  East  contiguos 
cell  has  index  1 ■ 1. 


9  All  cells  on  band  have  been  processed. 

(1)  ITA(l)  - ITA(2) 

(ii)  1X1(1)  - 1X1(2) 

(ill)  WXN(l)  - WNJ  - WXN(2) 

(iv)  WXS(l)  - WSJ  - WXS(2) 

(v)  WXN(2)  - P4(3,K1,I) 

(vl)  WXS(2)  - P4(3,K2,I) 

(vii)  DAI  - P4(1,K2,I) 

(vlli)  DA2  - P4(2.K2.I) (GTISE.  107) 


That  is  (1)  ^1,1,3  “ '’2,i-l,j 

(ii)  “ 52,1-1,3 
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1 


J 


(KMI) 


I 


where  is  one  of  the  3-D  prognostic  variables.  Note  replacement  of  data 
In  P2  which  Is  no  longer  needed. 


P2(2.K2.I)  - DA2 


(GTISE.127] 


All  grid  values  remain  the  same  for  ts.  N as  for  ts.  N - 1. 


IM  - I 
I - I + 1 


(GTISF.129] 


Process  next  eastward  cell. 


g (IPX.  EQ.  0)  GO  TO  1010  $ DO  1001  N - 1.  NBJ 


If  this  Is  not  a display  step,  return  to  ROTATE.  If  It  Is  a display 
step,  place  updated,  band  j data.  In  compressed  format  Into  the  array  BLK. 
Then  check  to  see  If  current  band  Index  J,  Is  the  last  band  of  a compressed 
data  block.  If  It  Is  not,  return  to  ROTATE.  If  It  Is;  then  buffer  BLK  onto 
display  file,  reinitialize  the  LT  Index  for  BLK,  and  Increment  the  counter 
for  the  number  of  I/O  blocks  processed  (LB).  Then  return  to  ROTATE. 

15  CALL  BERROR ( J . ICE . lUN) 

ICE  Identifies  the  routine,  J gives  the  current  Index  when  the  error 
occurred. 

Note;  We  only  present  the  EULER  routines  STIS<^,  TIS<^.  Other  Integration 
routines  will  differ  In  detail,  but  not  In  substance  from  that  presented 
here. 


SUBIODTINE  RHS(PJ.PJM.PJP.PHJ.PHM.PHP) 

This  subroutine  returns  for  column  (1,J)  the  right  hand  side  of  each 
of  the  prognostic  equations.  The  tendencies  are  passed  through  CEF  common. 
Coimson  Blocks 
/BUF/ 


/UNITS/ 


/FG/ 

/CEF/ 

/CEH/ 

/CEFG/ 

/FM/  Holds  forcing  terms  for  KV  levels  Indexed  (i,J). 

Variable  Description 

CORA  Array  returns  condensation  rate  for  all  KV  levels  from  CON  entry  of 
forcing  function  subprogram.  Presently,  values  returned  are  zero. 

CPRT  P*o/Cp 

ELPP  Array  of  L;^(P*(>)  at  the  KVl  level  boundarlett  (KH). 

ETPP  Array  of  L0(P*(fr)  at  the  KVl  level  boundarlea. 

FID  Array  returns  (F;^,Fg,F^,F  ) for  all  KV  levels  from  forcing  function 
subprogram  (KH).  Presently,  values  returned  are  zero. 

CLP  Gj^(F*) 

GLPT  Array  of  Gx(P*T)  at  all  KV  levels  of  current  index  (i,j). 

CTP  GgCP*) 

GTPT  Gg(P^T) 

H Array  holds  Hj^(Y)  or  H(K)  ■ 

IT  Index  for  n,  C arrays. 

IT  ■ 1 west  side 

IT  ■ 2 east  side 

KW  KV.  Avoids  control  of  do  loop  in  common. 

M Index  of  3-D  prognostic  variable. 

NU  Index  of  do  loop  for  north  v sum,  and  south  v sum. 

HUL  Lower  zonal  index  for  calculation  of  line  integral  along  v contiguous 

cells. 

HUNT  Code  set  if  calculation  proceeds  thru  maximum  cell  index  In 

this  case,  calculation  is  partitioned. 

NUNY  ■ 1 Proceed  in  2 steps 

NUNY  ■ 2 Proceed  in  I step 


k 

Z Hj^(l)^o  (KH) 


I 


S6 


1 


NDF  Upper  zonal  index  for  calculation  of  line  integral  along  v contiguous 
cells.  

NUSY  Code  set  if  calculation  proceeds  thru  maximum  cell  index  . In 
this,  case,  calculation  is  partitioned. 

NUSY  • 1 Proceed  in  2 steps 

NUSY  ■ 2 Proceed  in  1 step 

PHJ  ^ field  for  current  band  index  j. 

PHM  field  for  band  j - 1. 


PHP 

field  for  band  j + 1. 

PJ 

3-D  data  array  for  current  band  index  j . 

PJM 

3-D  data  array  for  band  j - 1. 

PJP 

3-D  data  array  for  band  j + 1. 

PRT 

P*o 

QDT 

Array  returns  added  heat  for  all  KV  levels  from  HEAT  entry 
function  subprogram.  Presently,  values  returned  are  zero. 

of  forcing 

T1 

P*o(f  + muj^) 

T2 

T3 

WT 

Temporary  variable  to  hold  W.,N 

^l»i»j  l«i»j  2,i,j 

Z 

Array  holding  vertical  P velocity  oi  at  all  25  levels  of  current  column 

ZBl 

u is  only  needed  diagnostically,  and  is  not  held  outside  RHS. 

0®) 

2(0^^  - 0k_^) 

ZB2 

2(0k+4j  - 0^.1,) 

I 

I 
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Helpful  Conments 


MB  - 1 


(GRHS.43] 


Line  integral  around  cell  begins  with  east-west  boundary.  These  are 
the  slatplest  to  process  since  * Wg. 

GLP  ■ .5*WEW  * (PH(l.Kl.IP)  - Pjg.Kl.JM))  $ L » IP 


P*E  ■ P*W 


where  is  the  value  of  P*  in  the  east  contiguous  cell,  (I  - IP)  and  P*y 
is  the  value  of  P^  in  the  west  contiguous  cell  (I  ■ IM) . The  contiguous 
cell  index  (L) , is  then  set  to  IP  in  anticipation  of  the  calls  to  YSL,  PST, 
and  PSP.  L is  passed  to  these  routines  in  FG  common. 

CALL  YSL(PJ.PT) 

Returns  Y1  array,  holding  | ^ — ~ all  KV  levels,  and 

YP  array  holding  or  all  KV  levels.  The  formal 
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arguments  PJ,  PJ  are  the  same  since  the  contiguous  cell  Is  on  the  same  band 
as  the  current  cell. 


CALL  PST(PJ.PJ) 


Return  PT  array  holding 
CALL  PSPCPJ.PJ.PHJ.PHJ) 


for  all  KV  levels. 


' if  ^ 

Returns  PH  array  holding  , 

l 2 J I Ok-y 


for  all  KVl  level 


boundaries. 


1 CONTINUE 


fGRHS.59J 


On  emerging  from  the  k indexed  loop 


Gx(P*T) 


To  -H  Tt 
2 


2 2 


D3(M,K)  - 


P*0  + Up  'i'o 


ELPP (KVl)  « PH(KV11 


L - IM 


^X^^*'^^ltvl 

(GRHS.63) 


(GRHS.611 
P*o  + P*E 


Integration  proceeds  along  the  western  boundary  of  a box.  Now  YSL  re- 
turns the  array  Y1  containing  at  all  KV  levels,  and 

the  array  YP  containing  ^*0  -0  '^'o  1"  at  all  KV  levels,  and 


PSP  returns 


f ^ 

the  PH  array  holding  (tp„  “ ) for  all  ^1 

I 2 °tJ5 


level  boundaries. 


3 CONTINUE 


(GRHS.78) 


On  emerging  from  the  k indexed  loop 


= l(P*«k  - “e»  I 


l±E  'Tq  Tr  _ P*o  P*V  Tq  + JW 


Lx(M)k  - ^ 


'k-lj 


^*o  P*W 
2 


Note:  the  factor  4 is  absent  from  (2*9  KH) . 
H(K) 


^ (1*3-11*^  . ^*o  ^*W  ^o  + -^W 

J 1 2 2 2 2 


D3(M,K)  - (!*£L±Is£ 


V 


2 j 


^*o  ^*W  f^Q  + 


ELPP(KVl)  - WEWJ  * {ELPP(KVl)  - PH(KVl)} 


+ ^F. 

2 

’<'0  »W  1 

2 / 

(GRHS.80) 

- ^ ) 
kv+Js  °kv+ii 


P*o  + P*5? 


"kv+ij  °kv+i5 


>} 


MB  - 2 (GRHS.88) 

Integration  proceeds  along  the  north  boundary.  The  easternmost  (L  ■ n^) 

contiguous  cells  are  handled  first.  Notice  in  the  calls  to  YSL,  PST,  and 

PSP,  that  PJM  identifies  the  northern  contiguous  band  on  which  the  adjacent 

cells  indexed  L are  located.  For  example  YSL  returns  ^ 

2 2 


in  Yl,  and 
loop , and , 


^★o  + p*m  Vo  + vr,. 

^ ^ ^ in  YP  on  the  first  pass  through  the 

2 2 2 

P*o  + P*n2.  ^ Yl,  and  ^*^2  '^>^1 

2 2 2 2 2 


in  YP  on  the  second  pass. 

5 CONTINUE (GRES.  114) 

On  emerging  from  the  loop 

W\  ■ ^24^ 
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8  Corresponds  to  case  v > 0 (KMl)  or  v„  ■ 1.  In  this  case  the  west- 
ernmost and  easternmost  cells  are  adjacent,  so  that  the  northern 
integration  Is  already  complete. 


9  Corresponds  to  case  v^-lorv^*!.  The  westernmost  and  eastern- 
most cells  coincide.  Its  necessary  then  to  subtract  the  value  above 
where  this  cell  is  counted  twice,  and,  then,  northern  Integration 
Is  complete.  Notice  in  this  case  that 

10  In  this  case,  there  are  v • n2  ~ (hj^+1)  cells  between  the  cell 
Indexed  nx  and  che  cell  Indexed  n2-  These  cells  are  indexed 
+ 1 thru  n2  “ I’- 
ll In  this  case,  Hx  “ ^i-x  *****  *^2  ^ **1-1*  **■*  v ■ ^ - 1 cells 

between  the  cell  Indexed  nx  and  the  cell  Indexed  n2-  These  cells 
are  Indexed  1 thru  r\2  ~ 


12  In  this  case,  no  • 1 and  nx  ^ *’j-l*  There  are  v ■ bj.x  “ cells 
between  the  cell  Indexed  nx  and'^the  cell  Indexed  n2-  These  cells 
are  Indexed  nx  + 1 thru 

13  In  this  case  nx  ^ ^’l-l  n2  ^ 1-  This  integration  is  handled  is 
two  stages.  First  TNUNY*!)  the  v - bj_x  ~ hx  cells  between  the 
cell  indexed  nx  and  the  cell  Indexed  1 are  processed.  These  cells 
are  indexed  nx  + 1 thru  bj«x-  Then  (NUNY*2)  the  v • n2  ” 1 cells 
between  the  cell  indexed  dj>x  and  the  cell  Indexed  n2  are  processed. 
These  cells  are  indexed  1 thru  n2  ~ 
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r 


8 CONTINUE 


(GRHS.1711 


On  emerging  from  the  north  NU  sum  we  have 

G8(P*)  - 

r /v  a.^  - ^ '^'>1  « ^ ^*°  ^*^2  "^o  + Tr,2 

2 2 *^'’1  2 ^"2 


+ E wN  1*0  - 
vjj  I 2 

^*o  ^*rn 

H ^ 


P*o  + P*!  To  + T, 


ni  ” '^o^k-is 


^*o  ^*rio 

* i ^ ‘‘"j  - ♦o>Vc-%  ^2 

+ z nKliS^iil.  (♦^  - , 

\K,  I 2 I o k-J* 


H(K)  - W™  /I^O-LZil.  - P*o  P*w  Op  ■*•  uw  1 

'^  1 2 2 2 2 1 


E*o  + P*m  Vo  + V|. 


Wni  + 


?*o  + E*no  Vo  + Vr 


+ E tjN  T*o  + P*\  Vq  ■*•  Vl 

Vj,  I 2 2 


I.K)  - P 


*o  + P*E  Vo  + ug  ^*o  •*•  'l*E  ^*o  P*w  Uo  + u;j  'l>o  + 


+ ^*o  ^*ni  Vq  •*“  vm  <>0  + w 


^»o  ^*n?  VQ  '**  V n?  ’^o  '*'  w 


+ Z 


P*0  + P*T  Vq  •*•  Vi  »o  •*•  »1 
2 2 2 


C SOOTH  BOUNDARY 


(GRHS.I73: 


The  south  boundary  Integration  proceeds  exactly  as  the  north.  Note  in 


particular,  that  the  contiguous  band  in  this  case  is  PJP. 
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n 


24  CONTINUE 


(GRHS.256) 


On  emerging  from  Che  south  NU  sum  we  have 

Oeff.)  - p.,  - 


- E W?  P 
''S 


1 

Vc  '■ 


gq(p  T).  hUliLu  +!i2JL!in2. 

0*k  2 2^1  2 2 '’2 


. P*o  + P*C2  '"o  ^52  '<'0  H2  c ^*0  + P*\  ^o  + V Ij/Q  + ip 

_ • • ^ - - 2 T-T  _ r T.tS  __  _ ^ 


2 


PRT  - PJ(1.  Kl.  I) 


"52  -v^s"? 


(GSHS.258) 


TRT  P^ 


-GTP  - .5  * (GTP  - WNSJ  * PRT)  (GRHS.260) 


- “»  «nV.nl-+  Jh"?'*!-  "u'’*u  - “u'Ha  - Ji'f  *, 


- f*„<%  - “s» 


ib:Qiii.2.8.  KH 


Ge(p*)  • (2^-  Eg)  Wi  - P*o<%  - %) 

- % - r ) {P*  W^  } + »s  p (Z^  - E3)  - P*  (W^  - Wg) 

- \ 


" ' f%  - Ws  Y 

- % (In  - + ^*o\;7-2 

- Ij  (Z^  - Zg)  - P*o^%  - above. 


GTPT(K)  « GTPT(K)  - WNSJ  * PRT  * PJ(3,k,I] 


CGRHS.263) 


G (P*T)^  - (I  - Z)  ^*°  ^ ^*1  '^o  ^ '^l  - (Z^  - Zg)  P*o 


GRHS.266 


H (1)  - H(l)  * SIGl  (1) 
Hid)  - Hi(l)  Ako 
IF  (KV  .EQ.  1)  GO  TO  43 


Skips  loop  If  only  one  level. 

H(K)  - H(K-l)  ■^  H(K)  * SIGl(K) 
k 

a,(K)  - Z Hi(K')  Aka  (3.6A  KH) 
kM 


RAP  • -H(KV) 


(GRHS.271) 

kv 

Z Hi  Aj^'o 
k'-  1 


(3.5A  KH) 


IF  (KV  .EQ.  1)  GO  TO  52 


(GRHS.273) 


If  KV  ■ 1,  then  u array  Is  complete,  and  equal  zero  at  top  and  bottom 
of  atmosphere. 
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2B(K)  - -RAP  * SIG(K)  - H(K-1 


Notice  that  ZB(K)  - 
kv 


f.o  “kj,  - Vs  ‘k'"  - ‘k'”  (3-6*  ™> 


5 * (ZB(K+1)  + ZB(K))  + SIG3(K)  * (RAIH<5LP*PJ(1.K.I))  .+  GTP  * PJ(2.K.J 


“k  “ ^*o 


IF  (KV  - 2)  53.  47.  48 


+ ^oGeCP*)]^ 


(3.7A  KH) 


Since  - (i^  - 0.  Z22^  - SIG2(KV)  * ZB(KVl)  - 0,  and  ZBlj^  - SIG2(1)  * 

ZB(1)  *0.  If  KV  ■ 1,  then  ZBl^^  ■ ZB2j^  • 0,  so  that  the  extra  term  (2.4  KH) 
does  not  appear  in  the  divergence. 

If  KV  • 2,  then  KV  - Jj  >j;  although  in  this  case-  then  ZBl^rr  ^ 0 and 


ZB2]^  t*  0;  so  that 


where 


M(M,KV)  - (Z 
E 


fPeo  + 

(Z  - Z)  ^ r 

E W I 2 


Peo  + P*\  uo  + ♦o  + ♦ 

2 2 2 


P»l  Vq 

2 2 


^*Q  '**kv-^  ’**kv  "^kv-l 

®kv+>i  “ ®kv-Js  2 


^kv-»i  " 


'^kv  ’l^kv-l 
2 


(2.4  KH) 


where 


■)3(M.i)  ■ (I  - 1)  » I 

E W I 2 2 2 

^ SM  2 2 2 


+ iio_^Z2_  ^ »2 

«3/2  - ^1/2  2 


"^3/2  - 


♦ l + <»2 


2 


95 


If  KV  > 2,  Chen  both  terms  ZBl  and  ZB2  are  present  simultaneously,  for 


indices  K » 2,...,KVM,  and 


D3(M,K)  - (I  - E) 

E W 1. 


*o  + uo  + '^o  *^1 


sM  2 2 2 


2 

+ Pi,  + „,  I'n  + 4», 


“>} 


*o 


■/_  '<'fc4-l  + 'f'k  - 

l“lc+is  ^ “k-ls  2 / 


CALL  FRC 
CALL  HEAT 


(GRES. 304) 


Forcing  function  subprogram 

RTE(l.K)  - -D3(1.K)  + T1  * PJ(2,K,I)  - R * GLPT(K)  - T2SIG2(K)  + FID(1,K) 


SUBROUTINE  YSL(PJ.  PJL) 

I This  routine  returns  Y at  the  coincident  boundary  between  cells  L and 
Jr.  at:  all  KV  levels.  The  version  here  described  Is  Che  filtered  factor  form; 
ocher  versions  of  YSL  have  been  programmed,  and  are  available  through  distinct 
update  deck  names  (GYSL2,  GYSL3). 

Common  Blocks 


/BUF/ 

/UNITS/ 


/CEFG/ 


Variable  Description 


Level  Index  for  do  loop. 

Index  for  3-D  prognostic  variable  In  YP  array. 


P*o  Pi 


P*0  + P*l  Yo  + 


where  Y is  u(MB  ■ 1)  or  v(MB  ■ 2). 


Helpful  Comments 


CCTITINUE 


(GYSL.391 


YSL  returns  the  arrays  Y1  and  YP  to  RHS.  For  Integration  along  the 


F*r>  + F*1  Ur,  + Ui 

east  or  west  boundary  (MB  • 1) , Y1  holds  — ° ^ ^ ° ^ — * 


and  YP  holds 


at  all  KV  levels,  where  i|i  Is  a 3-D  variable.  For 

2 2 2 

P 4*  P 

Integration  along  the  north  or  south  boundary  (MB  ■ 2),  T1  holds  ^ 


v„  + V, 


and  YP  holds 


P*o  + P*i  Vo  + V.  »o  + l<>i 


at  all  KV  levels. 
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SUBROOTINE  PST(PJ.PJL) 

:i"  This  routine  returns  PstarT  at  the  coincident  boundary  between  cells  L 
and ■ I , at  all  levels.  Other  versions  of  PST  have  distinct  deck  names 
;(GPST2),  the  version  here  described  Is  the  filtered  factor  form. 

Common  Blocks 
/FG/ 

Tiw 

TWifsT  “ ~ 

/AC/ 

/CEFG/ 

Variable  Description 
K Level  Index  for  do  loop. 


PRT 


P*o  + P* 


1 


Helpful  Comments 
1 CONTINUE (GPST1.20) 

PST  returns  array  PT  to  RHS,  For  Integration  along  a cell  boundary  PT 
P + P T + T 

holds  —^2 — 2—1 — L ^ where  I Is  the  zonal  Index  of  the  contiguous  cell. 


SUBROOTINE  PSP(PJ.PJL.PHJ.PHL) 

This  routine  returns  PstarPHI,  at  the  coincident  boundary  between  cells 
L and  I,  at  all  KV  vertical  levels.  The  version  here  described  Is  the  R and 
H fora,  other  versions  are  avallabTe  through  distinct  update  deck  names. 
Common  Blocks 
/FG/ 

/BUF/ 

/UNITS/  ~ 
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/AC/ 

/CEFG/ 


Variable  Description 


R 


PRT 


Level  index  for  do  loop. 


Helpful  Comments 
1 COHTINUE (GPSP1.37) 

P '4*  P 

PSP  returns  the  array  PH  to  RHS.  PH  holds  *°  — ^ (ij),  , ~ , ) 

2 ^k-»i  k-H 

at  all  KV  levels,  where  l is  the  zonal  index  of  the  contiguous  cell  on  its 
band  (JMl  for  north,  JPl  for  south,  J for  east-west). 


SUBROiniNE  ENTIS 

End  time  integration  scheme.  This  routine  completes  necessary  I/O 
at  the  end  of  a time  step.  INPUT  file  INI  and  IN2  are  rewound,  and  output 
file  lOT  has  records  JBl  to  JB5  written  on  it.  lOT  is  then  also  rewound. 
Common  Blocks 
/LTC/ 

/BUP/ 

/UNITS/ 

/AC/ 

Description  of  Variables 

1 Index  for  do  loop. 

NBJl  Number  of  cells  on  last  data  band  JB. 

Helpful  Comments 

2 REWIND  INI  $ NBJl  ■ NB(JBl)  $ CALL  CBUrT(IOT.0.99.U.PJI. I2*NBJ1) 

ENTIS  is  called  from  ROTATE  directly  after  the  call  to  TIS<^  was  pro- 
cessed with  J ■ JBl.  During  this  call  TIS^  updated  the  last  data  band  using 

W 


n 


P3(J*JB),  P4(J"JB1),  and  P5(J"JBl+l).  This  data  was  placed  In  array  P2, 
which  Is  passed  by  ROTATE  to  ENTIS  as  the  array  PJl.  Notice,  also,  the 
anticipatory  buffering  of  P6(J*JB1+2“JB3)  at  the  beginning  of  TIS(^.  For 
EULER  steps  then  the  south  polar  dummy  bands  JB2,  JB3  are  required.  For 
more  general  schemes,  dummy  bands  JB2-JB5  are  needed.  The  call  to  CBUFT 
writes  PJl  on  the  output  file.  For  the  single  level  case,  PJ  is  dimensioned 
4 X 3 X NBJX.  The  updated  JBl  data  is  in  4 x 3 x NBJl  > 12  x NBJl.  The 
99  is  a unique  identifier  for  the  BERROR  routine,  and  the  zero  identifies 
the  call  within  ENTIS. 

DO  3 I - 1,3 

3 CALL  PBUFT(I0T.I.99.  U.  BLK.3)  $ CALL  CBUFT(IOT.4.99.U.BLK..3) 

South  polar  dummy  bands  JB2-JB5  are  output  on  file  lOT.  PBUFT  performs 
a unit  check  before  the  buffer,  but  not  after,  and  thus  avoids  two  Initializa- 
tion times.  I,  again,  identifies  the  call  to  PBUFT  within  the  loop. 

IF  (IPX  .EQ.  1)  END  FILE  IDES  $ RETURN 

All  active  data  files  are  now  rewound.  If  a display  file  was  created 
in  this  ts.,  an  EOF  is  written.  AU  is  now  ready  for  generating  the  next  ts. 

APPENDIX  C ; DISPLAY 

PROGRAM  DISPLY 

Displays  files  saved  from  GEX  run. 

Characteristics 

Core:  125000^ 

Time : Open 

Options:  Described  below 


'1 

I 

( 


100  . 


/BUF/ 

101 

14600 

/UNITS/ 

14701 

15 

/AC/ 

14716 

2260 

Display  j 

^ REPORT 

17176 

657 

Options  { 

1 

DISPLY 

31451  . 

25564 

INITIALIZATION  1 

f SKFL 

57235 

220 

PROGRAMS  1 

[ GEN 

57455 

35 

SYSTEM 

f 

1 

102737 


12000 


Common  Blocks 


Parameters  used  In  spline  option. 

Axis  labels  and  parameters  for  plots. 
Used  to  expand  a display  file. 


/SPLNCO/ 

/PLC/ 

/BUF/ 

/UNITS/ 

/AC/ 

/DPC/  Parameters  used  in  spectral  plotting. 

Description  of  Variables 

AM  Array  of  second  derivatives  (moments)  returned  by  spline  fitting 

subroutine . 

AZLAB  Array  of  BCD  characters  used  to  label  x axis. 


AYLAB 

BLD 


CZ 


Array  of  BCD  characters  used  to  label  y axis. 

Array  of  values  of  current  variable  at  grid  points.  BLD  is  a mirror 
image  of  FLD,  since  the  CONTOR  subroutine  expects  Y values  in  in- 
creasing order. 

Array  of  spectral  components  returned  by  FFT  routine. 


Coefficient  of  constant  term  in  spline. 
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I 


Cl 


FILES 


INPUT 

(5) 

Community  Input  File 

OUTPUT 

(6) 

Community  Output  File 

INIT 

(7) 

Private  Binary  File  of  Grid  Information 

ISTARl 

(8) 

Scratch  File 

IDES 

(9) 

Display  File 

Display  Deck  (1  possibility) 


lTTACH,TAPE9,DX 

TTACh,  tape?, IX 


^ATTACH7 

^ NAME,  CM,  T. 


.Report  ts.  3360 
-Parameters  for  this  run 


.Compile  update  decks 


.Update  directives 


■Execute 


•Library  Card 
■Compile 


-Organize  Display  Decks 
-Attach  Plot  Routines 
-Attach  Display  File 
-Attach  Init  File 
-Attach  Global  Library 


C2  Coefficient  of  M term. 

C3  Coefficient  of  term. 

C4  Coefficient  of  term. 

DAYS  Simulated  time  In  days. 

F Value  of  specified  variable  around  band. 

FAC  Amplification  factor  In  plotting  run.  Now  1.0. 

FLD  Array  of  values  of  current  variable  at  grid  points. 

FLN  ■ hi  of  spline. 

NBR  Horizontal  and  vertical  wave  number  table  for  labeling  plot  of  2-D 

spectral  conponents. 

FTI  Simulated  time  Increment  At. 

H Normalization  factor. 

HD  Contour  level  Increment. 

HF  First  contour  level. 

HH  Last  contour  level. 

HL  Vertical  width  of  plotting  surface. 

lA  Index  used  to  organize  data  on  output  file. 

ID  30-character  plot  Identifier. 

lEND  lEND  * 0 Process  another  display  option 
lEND  " 1 Begin  end  processing 


IFLD  Variable  code  IFLD  Variable 

1 u 

2 V 

3 T 

4 r 

5 P* 

6 


IN  Location  of  current  variable  within  variable  field  on  IDES  file. 

INC  Increment  to  next  appearance  of  current  variable  on  IDES  file. 

INTRP  INTRP  ■ 0 Spline  option  not  called 

INTRF  ■ 1 Spline  routines  are  called 

1D4 


1 


} 

I 

I 


I 

I 


IPLT 

ISC 

11 
IIP 

12 
12? 
JB 
JP 
J1 
J2 
J3 
J4 
L 

LEV 

LT 

LTS 

LTSL 

LZ 

M 

NBJ 

NBJM 

NCL 

NLCX 

NLCY 

NP 

NPI 


IPLT  * 0 Run  does  no  plotting 

IPLT  ■ 1 Plotting  subroutines  are  called 

Scratch  file  name  (not  used) . 

Counter  used  to  organize  data  for  dumps. 

11  + 1 

Counter  uaed  to  organize  data  for  dumps. 

12  + 1 

Total  number  of  latitude  bands. 

Wave  number  J - 1. 


Loop  controls. 


Array  Index. 

\ 

Vertical  level  number  of  current  variable. 

Index  for  stripping  BLK  of  current  variable 
Current  time  step  displayed. 

Previous  time  step  displayed,  or  -1  If  no  previous  step. 
Number  of  words  In  a data  block  In  IDES  file. 

Loop  Index. 


Number  of  contour  lines  In  a plot. 

Number  of  BCD  characters  In  label  for  horizontal  axlr. 
Number  of  BCD  characters  In  label  for  vertical  axis. 
Number  of  pages  to  accomodate  data  dump. 

Loop  Index. 
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I 

I 


i 

I 

[■ 

I 


f 

II 

i 


NPX 

_Lpop  indi 

NT  T -- 

.2  * JB 

ST7- 

. Defines  ] 

NTY 

C*I 

NTY 

NTY 

NTY 

Cc.2T  r t 

NTY 

NTY 

* £ • £T.  « 

NTY 

HDHB 

_^ray  of 

N1  - 

N2  + 1 

N2 

JsNBJX 

Report 

Raw  data  dump  of  field  IFLD 
Interpolated  data  dump  of  field 
Plot  field  IFLD 
.Spectral  dump  of  field  IFLD 
Normalized  spectral  dump  of  field 
Spectral  plot  of  field  IFLD 


N2P  N1 

N3  Loop  index. 

RFLN  AXj 

T Output  string  from  contour  plot. 

X Zonal  grid  table. 

XL  Maximum  length  of  plotting  surface. 

XNP  Fractional  number  of  pages  for  dump  output. 

XNT  Truncated  number  of  output  pages. 

XT  Temporary  array  of  interpolated  values. 

7 Meridional  table  in  degrees. 

YT  Temporary  array  of  Interpolated  values. 

ZMAX  Maximum  of  current  variable  over  globe. 

ZMIN  Minimum  of  current  variable  over  globe. 

Helpful  Comments 
ISTARl  - ISC 

ISTARl  ■ 8 is  used  as  a scratch  file  to  expand  compressed  data  file  for 


'1 


ll 


i 

) 

'i 


report  output 
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CALL  CBUFN(INIT.0.0.U.CBLK.1200) 


(GSPLAY.25) 


Fills  out  AC  common  with  grid  Information  for  display  file. 

CALL  GEN(X.Y)  $ IF  (INTRP.  EQ.  1)  CALL  SPREP(NBJX) 

GEN  returns  arrays  X,  Y;  ■ (l-l)^^AXj } and  y^  ■ ®jb1-J 

J 

degrees.  SPREP  Intltallzes  values  for  the  cubic  splines  routine.  NBJX  has 
no  slgnlflcwce  in  the  call  here;  and  Is  not  used  in  SPREP. 

AXlJbrY)  « SH_DAYS fGSPLAY.30) 

Axis  labels  for  contour  plots. 

-I-  MP  * (Kl-1) 

Determines  number  of  variable  position  to  pass  over  on  compressed  data 

file. 

LTS  - LTSL 

Initializes  last  time  step  accessed  on  display  file  to  -1. 

699  IF  (LTS.  GE.  0)  LTSL  - LTS 

If  a display  file  vas  previously  accessed,  save  Its  ts.  In  LTSL.  This 
value  will  be  used  to  seek  the  next  file  (LTS)  in  the  most  economical  way. 

All  display  options  will  loop  back  to  699  for  a decision  on  further  processing. 
IF  (TEND.  EQ.  0)  GO  TO  6999  $ IF  (IPLT  .EQ.  1)  CALL  ENDPLT  $ STOP 

Processing  ends  If  the  previous  read  statement  set  the  code  lEND  ■ 1. 

A call  to  ENDPLT  Is  mandatory  If  this  run  created  any  plot  files. 

6999  READ (INPUT. 700)  LTS.  IFLD.  LEV.NTY.  TEND 

Time  step  LTS  Is  to  be  accessed  on  the  IDES  file  and  displayed.  The 
nature  of  the  display  is  determined  by  the  values  read.  See  sample  deck. 

DAYS  - LTS  * FTI/ 86400. 

Simulated  time  of  ts.  LTS  In  days. 

ENCODE  (10,502)  AXLABO))  DAYS  $ CALL  SKFL (IDES. LTS. LTSL) 


These  characterizing  numbers  are  refotwtted  In  internal  BCD  code,  to  be 

-TOT 


I 


used  as  axis  labels  if  plotting  Is  Indicated.  The  call  to  SKFL  locates  the 
display  file  labeled  LTS,  and  positions  that  file  directly  after  the  3-word 
label. 

IF  ((NTY  .NE.  0)  .OR.  (IFLD  .NE.  0))GOTO  698  $ BACKSPACE  IDES 

The  code  for  report  output  has  both  NTY  ■ IFLD  ■ 0.  In  this  case  the 
IDES  file  LTS  must  be  first  expanded  onto  a scratch  file  (ISTARl) , before 
It  Is  sent  to  REPORT.  Since  EXPAND  expects  a file  positioned  at  the  beginning 
of  the  3-word  label,  the  backspace  over  1 record  Is  done. 

GO  TO  (1000.  2000.  3000.  4000.  5000.  6000)  IFLD 


The  axis 

labels  are 

now  ready 

1 

2 

3 

4 

AXLAB 

I 

T 1 

OZ  > 

ts  # 

days  # 

DAYS 

AYLAB 

DATA  SER 

serial  9 

GLOBAL  RU 

N 

The  value  of  IN  locates  a variable  within  a data  block  on  a compressed 
file.  P^  and  occupy  positions  1 and  2 respectively,  then  u,  v,  T,  r fol- 
low consecutively  up  a column.  IN  locates  one  of  these  variables  at  a fixed 

level. 


The  logic  of  the  loop  Is: 


N3  - 1,  NBLK 

LT  * IN  Locate  first  occurrence  of  variable  within  present 

data  block  at  level  to  be  ripped  off. 


J1  • J2  + 1 First  band  Is  present  block 
J2  > JLI(N3)  Last  band  Is  present  block 

LZ  ■ LSZ(N3)  Number  of  words  In  present  block 

, Buffer  In  data  block 

BUF  (LZf* jj3 

J - Jl,  J2 

NBJ  ■ NB(J)  Number  of  cells  In  present  band 
L • J - 1 Latitude  Index  in  FLD  array 

I - 1,  NBJ 

FLD(I,L)  • BLK(LT)  Rip  off  variable  within  present  data  block. 

LT  ■ LT  + INC  Next  occurrence  of  variable  at  sane  level  Is  separated 


by  INC  words. 


6002 
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-Notice  that  before  entering  the  loop  IDES  Is  positioned  after  the  3-word 


label  for  ts.  LIS.  On  leaving  the  loop,  IDES  is  positioned  before  the  EOF 
of  ts.  LIS,  and,. FLD  contains  a full  globe  of  the  specified  variable  at  a 
fixes  level.  In  the  form 

M 

iF~'(NTY  .NE.  0)  GO  TO  6003  $ XNP  - FLOAT(NBJX)/8.0  $NP-XNP$XNT»NP 

N^  « 0 produces  a raw  data  dump  of  the  field  held  In  FLD.  The  printed 
output  should  write  the  two  dimensional  array  with  a flexible  number  of  cells 
on  each  band.  The  data  Is  to  be  reported  horizontally  by  zonal  Index,  and, 
vertically  by  latitude  Index.  To  print  page  numbers,  the  data  Is  partitioned 
by  determining  the  number  of  pages  needed  to  report  the  band  with  th4  largest 

I 

number  of  cells,  eight  cells  to  a page  (XNP).  This  value  Is  truncated  (NF), 


FLD 


and  then  floated  (INP) . If  XNP  > XNT,  then  reporting  8 cells  to  a page  will 

require  an  extra  fraction  of  a page  to  accomodate  the  band  with  NBJX  cells 

(NP  - NP  + 1). 

The  logic  of  Che  loop  Is: 

II  - -7 

NPX  - 1,  NP 

11  * J1  + 8 Smallest  zonal  Index  on  this  page 

WRITE (IPRT, 7001)  Write  specified  variable  (labeling  and  num- 

bering pages). 

12  ■ II  -t*  7 Largest  zonal  Index  on  this  page 
J - 1,  JB 

NB  ■ NB(J+1)  Number  of  cells  on  data  band  J 

WRITE  NUMB ( JIP )— NUMB (I2P)  NUMB(M4-1)  - M so  this  writes  Che 

zonal  Indices  II  thru  12  across  page 
12  ■ MIN^(Il-f7,NBJ)  There  are  a maximum  of  NBJ(fNBJX)  cells  to 

be  reported  on  this  band 

^ n > 12? 

WRITE  NUMB(J+1)  Write  latitude  Index  J followed 

FL0(I1,J)  thru  FLD(I2,J)  by  segment  II  thru  12  of  variable 
GO  TO  6004 

> WRITE  (IPRT, 7003)  Skip  record  since  all  cells  on  this  band  have 

been  written 


004 
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The  situation  II  > 12  can  occur  i£  data  thru  NBJ  has  been  reported  on  page 


' 


NP  - 1 (so  that  12  - NBJ).  Then  Iljjp  - + 8,  so  that 

> ■ ..  -1  nbjx 

t I , - I . , I 

II  I NBJ  Ilup  11+7 

reported  on  page  NP  - 1 

then  12  ■ MlNcf (11+7, NBJ)  < Tljjp.  Indicating  data  band  j reported. 

6003  IF  (NTY  .EQ.  1)  GO  TO  6011 
IF  (NTY  .GT.  2)  GO  TO  6020 

NTY  ■ 2 indicates  a contour  plot  of  the  variable  in  FLO.  Such  plots 
can  then  be  made  level  by  level,  for  any  specified  variable.  The  array  FLO 
Is  first  flipped  into  BLO,  so  that  column  1 of  BLD  contains  data  next  to 
the  south  pole,  and  column  JB  contains  data  next  to  the  north  pole.  Since 

the  Y array  returned  by  GEN  is  in  increasing  order  (9g  ■ -90° 6jj  • +90°), 

the  variable  field  must  be  in  the  same  order.  That  is,  J ■ 1 should  corre- 
spond to  the  grid  point  next  to  the  south  pole,  and  J > JB  that  next  to  the 
north  pole,  as  in  array  BLD. 

CALL  CONTOR  CBLD.66.34.X.Y,NBJX.JB.XT.YT.T.998,HF.HD.HH.NCL) 

Produces  contour  plot  of  data  field  BLD. 

CALL  PLOT(13.0. .-3) 

Moves  pen,  with  point  up,  13  inches  to  right  of  origin. 

C INPERPOLATE  TO  REGULAR  GRID 

In  this  option,  the  data  for  the  variable  in  FLD  is  interpolated,  band 
by  band,  to  a regular  grid  of  NBJX  data  points  per  band.  Cubic  splines  are 
used  to  replace  the  discrete  data  around  a band,  with  continuous  polynomials, 
on  which  the  data  is  evaluated  at  the  NBJX  regular  argtiments. 


no 


-6014  F(I)  - FLD(I.J)  $ CALL  SPREP2(NBJ)  $ CALL  SPCOFP(NBJ) 

F holds  the  NBJ  grid  values  of  the  variable  around  the  current  band  J. 
SPREP2  evaluates  certain  constants  needed  for  the  matrix  elimination  for 
this  band,  and  SPCOFP  produces  the  coefficients  of  the  cubics  fitting  this 
'data,  along  with  the  second  derivatives  at  the  grid  points. 

CALL  SPINT(NBJX.DEV.FLD.NBJ.J) (GSPLAY.88) 

Replaces  column  of  FLD  with  NBJX  interpolated  data  values. 

X INTERPOLATED  DATA  DUMP (GSPLAY.90) 

The  Interpolated  data  is  reported  horizontally  by  zonal  index,  and 
vertically  by  latitude  index.  Pages  are  numbered  as  described  previously. 

There  is  no  need  for  the  test  II  > 12,  since  all  columns  of  FLD  are  of 
length  NBJX. 

C DOUBLE  PERIODIC  FOURIER  TRANSFORM (GSPLAY.98) 

Spectral  data  options  follow. 

6020  CALL  DPFT 

Returns  2-dimenslonal  spectral  components  in  FLD. 

6036  WRITE  IPRT,712)  JP.  (FLD(I. J) .I-I1P.I2P)  $ IF  (NTY  .EQ.  3)  GO  TO  699 
The  loop  above  this  point  points  the  2-dimensional  spectral  components 
of  the  field  FLD,  numbering  pages,  and  wave  number  components  horltontally 
and  vertically. 

6045  CONTINUZ  $ IF  (NTY  .EO.  4)  GO  TO  699 

The  loop  above  this  points  the  2-dimenslonal  normalized  spectral  com- 
ponents of  the  field  FLD,  numbering  pages  and  components. 

6060  NLCX  - 50  $ CALL  CONTOR(FLD,66,FNBR,FNBR,NyBi;cr,YT,T,998,HF,HD,HH,NCL)  $ 
NLCX  - 40  $ CALL  PLOT  (13.0. .-3)  $ GO  TO  699 

Plots  field  of  2-dlmensional  spectral  components  normalized  above. 


111. 


6050  BACKSPACE  IDES  $ REWIND  ISTARl  $ CALL  EXP AND (IDES, ISTARl) 

CALL  REPORT (ISC, -1)  $ GO  TO  699 

The  -1  codes  REPORT  to  just  print  the  energies  and  masses  for  the 

specified  time  step. 

SUBROUTINE  GCUX.Y) 

Returns  to  DISPLY  arrays  X,  Y.  These  are  the  zonal  and  meridional 
tables  passed  to  the  plotting  routines,  if  that  option  is  chosen  in  DISPLY. 
Common  Blocks 
/BUF/ 

/UNITS/ 

/AC/ 


Variable  Description 

DX  Minimal  zonal  Increment  in  degrees 
J Data  band  index. 

JJ  Reverse  index. 

RTD  ^ conversion  factor  from  radians  to  degrees 
Helpful  Comments 


. 180 
. If  . 


10  X(I)  - X(I-l)  + DX 

X(l)  - 0,  . . .,  X(K)  - (K-l)AXj^,  . . .,  X(NBJ)  - N(BJX-l)AXj^ 
Thus  X(NBJX+1)  i X(l)  - 360° 

11  Y (5)  - RTD  * (BAND(JJ^-3)  BAND(JJ-»-2)) 

Y(l)  - RTD  * (BAND(JB2)  + BAND(JBl))  - RTD  * Ojb+j,  + ®JB-»5^  " ®JB 


Y(JB)  - RTD  * (0,  + 0,  ) - 01 

■*'2  /2 

Then  Y(%)  5 0 , is  the  south  pole  (-90°)  and  Y(JB+^j)  = 0j^  is  the  north 

pole  (+90°).  Note  GEN  returns  X and  Y in  increasing  order  for  CONTOR. 


I 


SUBROUTINE  SKFL(IU.L.LTSL) 

Subroutine  SKFL  is  called  from  DISFLY  to  seek  the  file  labeled  with  the 
ta.  L,  on  the  tape  lU.  On  entry,  lU  is  positioned  either  at  the  beginning 
of  information  (LTSL  * ~1),  or,  at  the  end  of  the  data  block  of  the  time 
step  (LTSL)  previously  found  by  SKFL. 

Common  Blocks 

' ' ' rsfi 1 

/BUF/  — — - 
/UNITS/ 


Variable  Description 

N Do  loop  index  and  code  passed  to  CBUFN  routine  for  locating  source  of 
buffer  error. 

NBL  Value  sec  if  lU  has  to  be  backspaced  to  locate  L (LTS\<  LTSL) 

Helpful  Comments 


IF  (LTS  .LT.  0)  GO  TO  21  $ NBL  - 2 * (NBLK  + 1)  + 1 


GSPLAY.141 


If  LTSL  ^ -1,  then  lU  is  located  at  Che  end  of  Che  data  block  of  Che 
file  LTSL.  To  access  Che  label  of  Che  preceding  file,  we  must  skip  over  2 
data  blocks  (each  NBLK  records) , 2 labels  (each  a single  3 word  record) , and 
1 end  of  file  mark  (a  single  record).  Thus  NBL  is  set  to  2 * NBLK  + 2 -t*  1. 

1 BACKSPACE  lU  $ RETURN 

In  this  case  L « LTSL,  and  we  must  simply  backspace  over  the  data 

block  (NBLK  records)  to  put  lU  at  the  end  of  the  correct  label. 

IF  (ITS  .LE.  LD)  GO  TO  15  $ NBL  - NBLK  + 3 $ GO  TO  10 
15  WRITE  (IPRT.  100)  L.LD 


If  ITS  ^ LD,  Chen  L has  not  been  found  to  follow  the  first  file  on  lU 

(time  step  LD).  An  error  is  indicated.  Otherwise,  we  continue  backspacing 

over  time  steps,  and,  since  lU  is  located  at  Che  end  of  Che  label  for  ts. 

ITS,  to  access  the  label  (of  the  preceding  file  we  must  skip  over  1 data 
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block  (NBLK  records),  2 labels  (each  a single  3 word  record),  and  1 end  of 
file  mark  (a  single  record).  Thus  NBL  = NBLK  +2+1. 


, 

I 


IF  (ITS  .GE.  LDf)  GO  TO  40  $ DO  3 N =■  1,  NBLK 
3 CALL  CBUFN(IU.I.N.U.BLK.3)  $ GO  TO  20 (GSPLAY.153) 

If  ITS  ^ LIM,  then  L has  not  been  found  to  precede  the  last  file  on  lU 
(time  step  LIM).  An  error  is  indicated.  Otherwise,  one  data  block  (NBLK 
records)  is  skipped  over  to  prepare  to  access  the  label  of  the  file  following 
ITS  on  lU. 

SUBROUTINE  SPREP(N) 


This  routine  returns  to  DISPLY,  through  SPLNCO  common,  the  coefficients 
of  a cubic  spline  fit  to  the  data  in  variable  F,  and,  also,  the  second  deriv- 
atives at  the  fitted  points.  F holds  a full  band  of  NBJ  data  points  in 
SPLNCO  common.  The  procedure  is  described  in  ANW.  In  the  notation  found 
there,  the  spline  is 


Sa(x)  - Mi-1 


Mi 


(x  - x^_^)' 


6hi 


L , . 1 

f ^ 

“1  - * 

Mihi2 

'x  - - 1 ■ 

Fi-l ^ 

1 bl 

'T 

" 6 

S.  / 

6 

where,  is  the  second  derivative* 

Yi  - Fi 

hi  is  the  grid  spacing,  here  AAj 
For  the  periodic  case:  ■ M^j^ 


yi  • fn+i 


For  equal  spacing: 


jLt±l. 


AXj 


hi  + hi+i  hi  + hi+i  2AXj 


- H 


so  that  the  defining  equations  for  the  second  derivatives,  M^,  are: 


“Oi  X w 

[Mil 

pi 

u Di  X 

M2 

<12 

• 

• 

■ 

' U X 

• 

• 

_X  M Dl 

-Mn. 

L^J 
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The  first  step  is  to  eliminate  between  the  first  and  last  rows. 


A 

u 

X 

Ml 

• 

m 

1 

1 

□ - ^1^ 

• . . u(Di-PIQi)_ 

• 

• 

- Pldl_ 

J 

with  * U,  PIQj^  “ ^1^1*  eliminating  between  the 

first  and  second  equations: 


"Di  a Qi 

Mf 

di_ 

D2  A Q2  . 

• 

^2  ""  ^2^  1 

_ -Pj^A  . . . u(D2^-PIQ2^)_ 

- Pl<ll_ 

with  A2  * y/Dl.  D2  * - A2I,  Q2  “ "^2^1’  (below)  * -Pj^X/D2, 

PIQ2  - PIQi  - P2Q2- 

Then  eliminating  between  the  second  and  last: 


~Dl  A 

D2  A 

Ql  “ 

Mi" 

dl 

d2  **  A2dj_ 

-P2A 

uP1Q2_ 

_djj-P  2.*^!“^  2 ^d2“A2d2^)_ 

and  so  on  until: 
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X 

D2 


Qi 

^2 


Vz  ^ "^"”2 


-Ml' 

• 

r 

«o  *0 

1 

• 

a 

• 

• 

^N-1 

,1 

N-2 


where  % - D^  - E (PIQ)^,  5^  - (d^  - . .)  or,  finally, 

1^1  p _ ij 

eliminating  between  the  last  two  equations  using  ■ ^ - -p  , - A., 

Vi  ^ 


Di  X 


Dz 


®N-2 


Ql 

QZ 


%-2 


«r 

• 

1 

1 

• 

- 

h 

• 

• 

• 

•^N-l 

-1 

-MnJ 

uS  J 

where,  Ai  - w/D^.i,  - AiX,  ^i  - -Pi-iX/D^, 

I 

PIQI  - E PjjQj^,  then  + Q.Mj,  - 6^,  1 - 1 , 2,  . . . , N-1 

k"l 

and 


C.  S 

Common  Blocks 


/SPLNCO/ 

Variable  Description 

A Array  A(I)  ■ u/D(I-l)  used  as  multiplier  in  elimination  procedure. 
AIM  X - >s 

AM  Moments  M(I)  passed  in  SPLNCO. 

AMD  u • i] 

Cl  First  coefficient  of  cubic. 

C2  Second  coefficient  of  cubic. 


HI 


C3  Third  coefflclenc  of  cubic. 


C4 

D 

% 

F 

FAC 

FLN 

GI 

GIF 

I 

II 
IM 
IP 
N 

NM 

P 

PIQ 

Q 


Fourth  coefficient  of  cubic. 

Array  D(I)  of  final  diagonal  entries  of  matrix  after  elimination, 
diagonal  element  is  special. 

diagonal  element  in  last  row  after  elimination. 

Array  F(I)  of  data  around  band  to  which  spline  is  fitted. 

hi  ■ equal  intervals. 

AAj 

Xj^^(C4xj^  + C3) 

Index  for  loops. 

Index  accesses  arrays  from  bottom  to  top. 

I - 1 
I + 1 

NBJ  passed  by  DISPLY. 

NBJ  - 1 

Array  P(I)  of  multipliers  in  elimination. 

Array  PIQ (I). 

Array  Q(I)  of  entries  in  last  column  after  elimination  procedure. 


RFLN  AXj 

S Last  component  of  right  hand  sector  after  ellmlnatloa. 

TPI  2it 

XI  - (I  - l)AXj 

XIP  XI  + AAj 


Last 
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Heloful  Comments 


SUBROUTINE  SPREP(N) 


(GSPLN. 14) 


First  entry  calculates  the  arrays  A,  D,  Q,  P,  and  PIQ  for  the  elimination 
process.  Since  these  arrays  will  remain  the  same  whatever  the  right  hand 
side,  SPREP  need  only  be  called  once  at  the  beginning  of  DISPLY.  N is  not 


used  here. 


ENTRY (SPEP2) 


(GSPLN. 34; 


This  entry  calculates  values  which  appear  in  the  elimination  procedure. 

2 

In  particular,  the  value  of  N ■ NBJ  is  used  to  calculate  l/AXj,  3/AXj  and 
Dj|,  which  appear  in  equations  A. , B.  , and  C.  respectively. 


DN  - D(N)  + Q(N)  - P(NM)  * ALM  - PIQ(NM) 


Dn  ■ % + Qtl  - Pn-1^  ■ ^^^N-l 


(GSPLN, 40) 


ENTRY  SPCOFP 


(GSPLN. 42) 


This  entry  computes  the  coefficients  of  the  spline,  and  the  second 


derivatives. 


F(N+1)  - F(l) 


(GSPLN. 44) 


periodic  case. 

) - (F(2)  F(N)  - 2 * F(l))  * FAC 


(GSPLN. 45) 


“ II7  <^2  ^N  - 2yi)  - 

AM(I)  - (F(IP)  F(IM)  - 2.  * F(I))  * FAC  - AM(IM)  * A(I) 

Accumulates  the  values  of  the  right  hand  vector  in  the  AM  array.  That 
is  • d^  - in  particular: 

- dj,  - Ajj(dj,.j_  - A{,_2_(dj,_2-  • •)) 


S - S P(IM)  * AM(IM) 


(GSPLN. 49) 


Accumulates  the  values  of  the  last  component  of  the  right  hand  vector. 
That  is,  on  leaving  the  loop:  ~ 


i 


r 

j 


I 


i 

! 

I 

( 

1 

I 

I 

i 

f 


I 


i 


S ■ P + . . .;  the  actual  value  of  the  last 

component: 

® “ ‘^N  “ ^^‘^N-l  ” ^-l^‘^N-2  ■•••))"  ^l^^l  ■ ^2  ^**2  “ ^2*^1^'  • • 
■ (M^  - S)  above. 

AM(N)  ■ (AM(N)  - S)/DN (GSPLN.50) 


Mjj  ■ S/D^  gives  actual  value  of  Mj|. 


3 AM(I)  « (AM(I)  - ALM  * AM(IP)  - Q(I)  * AM(N))/D(I) 

Loop  solves  for  the  i ■ N-1,.  . .,1.  .We  already  have  Mjj. 

From  C.  » (<Sj^  - ~ we've  accumulated  5^  in 

X(I)  ■ XIP (GSPLN.58) 

Xi  - (i-l)dXj 

4 Cl(I)  - F(I)  - GI  - C^(I)  * X(I) (GSPLN.65) 


Since  from  A. 

SA(x) 


6AXj 


(xi+1  - x)^ 


Mj+l 


(x  - x^)^  + 


[Mi^^l  1 

'xi+i-x 

6 j 

. 4Aj 

then,  expanding  this  and  replacing  by  + AAj 

SA(x)  « C4(i)x3  + C3(i)x2  + C2(i)x  + ^(i) 

SUBROUTINE  SPINTCM.D.FLD.N. J) 

SPINT  returns  to  DISPLY  the  array  FLD,  which,  previous  to  the  call,  had 
only  NBJ  ■ N grid  values  in  its  column.  These  values  are  the  cubic 
spline  interpolated  values  around  the  data  band,  calculated  so  that  all 
bands  will  have  NBJX  interpolated  values  on  the  same  longitudinal  slices. 
Common  Blocks 
/SPLNCO/ 


12a 


Variable  Description 


ARG  X ■ OjAXj^  + (l-l)AXj^  where  Oj  is  the  fractional  deviation  of  the 
data  band  (Indexed  j+1). 

D Oj  passed  from  DISPL7. 

DARG  AAj^  ■ 2ir/NBJX,  where  is  the  index  of  the  band  with  maximum  cells 
NBJX. 

FLO  Data  field  to  be  interpolated. 

L Location  of  Interpolated  argument  x in  successively  Indexed  splines, 

if  the  splines  are  indexed  east  from  the  zero  meridian. 

M NBJX 

N NBJ  for  current  band. 

Helpful  Comments 

DARG  " TPI/FLOAT(M)  $ ARG  - D(J+1)  * DARG  $ L - 1 
- DARG  - 2it/NBJX 

The  first  argument  is  x ■ a.AA.  . This  places  x on  a reduced  scale  within 

J Jx 

its  own  first  displaced  box.  Since  AAj^  < AAj  this  x is  clearly  within  the 
first  Indexed  spline  L ■ 1. 


FLD(I.J)  ■ ((C4(L)  * ARG  C3(L))  * ARG  + C2(L)  * ARG  ClfL) 

Displaces  the  1^^  data  value  on  band  J with  the  interpolated  data  value. 


ARG  « ARG  4-  DARG  $L«FLN*ARG$L«  MOD(L.N)  ->•  1 


Next  argument  is  x + AAj^.  For  this  x:  L 


AA^ 


+ 1 gives  the 


location  of  the  argument  in  the  indexed  splines.  L • K occurs  if  J • Jx. 


Then  MOD (NBJX, NBJX)  + 1 - 1 is  ok. 

SUBROUTINE  CONTOR(Z.IZ.JZ.X.Y.M.N.XT.YT.T.K.HF.HD.HH.NTAPE) 


Produces  contour  plot  of  the  function  Z(X,Y),  which  Is  given  for  the 
grid  points  (X(I),Y(J)),  I - 1,.  . . ,M  and  J ■ 1,.  . . The  output  is  a 
sequence  of  strings,  for  the  plotter,  in  the  format  (L,H)  (X(l) ,Y(1) ,X(2) , 
Y(2),.  . .,X(LL,Y(LL))  where  L ■ 2 * LLapd  L K.  The  output  is  terminated 


by  a string  having  L ■ -1.  Each  string  describes  a contour  (or  part  of  a 
contour)  at  the  level  Z ■ H. 

Common  Blocks 
/PLC/ 

Variable  Description 

BLANK  Option  to  set  certain  z values  to  5HBLANK  so  that  they  may  be  skipped 
in  processing.  Not  used  here. 

DX  X increment  returned  by  plotting  routine  SCALE  (which  scales  x array 
to  fit  over  axis  of  10  Inches. 

DY  y increment  returned  by  plotting  routine  scale. 

FCT  An  amplification  factor  for  characters  LLINE  ^ 13. 

H Current  contour  level. 

HD  Increment  between  contour  levels. 

HDD  Increment  within  current  contour  level  ■ HD/ 1000 

HF  First  contour  level  ■ x^iin  *5-  ^*mln 

HH  Ust  contour  level  - < 5max 

HLEVEL  Current  level  or  preceding  level 

IFLAG  IFLAG  - 1 CONTINUE  SEARCH  FOR  LINKS 
IFLAG  - 2 OUTPUT  STRING  READY 

II  Determines  search  area  for  neighbors  in  link. 

IL  Zonal  index  of  preceding  link  point. 

ISYMB  Not  used. 

IS2  Index  for  loops. 

IX  Current  length  of  output  string. 

IZ  Zonal  dimension  of  z array. 

II  Indices  K1,K2,K3,K4  during  search  for  neighboring  points. 

147  147*1  Reached  mazlmum  size  of  output  string 

147*2  Continue  processing  ~ 

JJ  Determines  search  area  for  neighbors  in  link. 


JL  Meridional  Index  of  preceding  link  point. 

JLINE  JLINE  * 1 produces  a line  plot  with  symbols  at  every  point. 

Meridional  dimension  of  Z array. 

X”  Maximum  size  of  output  string. 

KEY  Code  passed  to  subroutine  NEIBOR 

KEY  ■ 1 Current  point  on  a vertical  link. 

KEY  ■ 2 Current  point  on  a horizontal  link. 

KIKI  Keeps  count  of  number  of  HDD  Increments  that  are  added  whenever 
current  level  passes  through  function  value. 

KL  Value  of  KEY  for  preceding  point. 

K1 

K2  ; • 

Codes  search  for  neighboring  points. 

K3 

K4  . 

LLINE  A number  describing  the  centered  symbol  to  be  used  1 £ LLINE  ^ 13. 

LSYMB  Array  holds  list  of  plotting  characters  used  for  table  to  right  of 
plot. 

M Length  of  X table. 

MMl  M - 1. 

N Length  of  Y table. 

NLABX  Passed  value  is  NLCX  « 40  from  DISPLY.  Determines  number  of  charac- 
ters in  X axis  label. 

NLABY  Passed  value  os  NLCY  > 40  from  DISPLY.  Determines  number  of  charac- 
ters in  Y axis  label. 

NMl  N - 1. 

NNLABX  HNLABX  * -NBLABX.  Sign  determines  orientation  of  characters,  clock- 
wise or  counterclockwise  to  axis. 

NSTEPS  Passed  value  NCL  from  DISPLY.  Number  of  contour  lines  in  plot. 

N1 

Determine  search  area. 

N2  . 
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r 


N3 


Determine  search  area. 


N4 

SKALE  Scale  factor  In  plotting  (now  1.0). 

SYSZ  Height  of  characters  plotted  by  SYMBOL. 

T Output  string  temporary  array. 

^ ^ " I’Wx  ^max  points  already  processed. 

X Table  of  zonal  grid  points  In  degrees. 

XF  Zonal  value  of  first  point  of  a string. 

XMIM  Minimum  X after  scaling  to  10  Inches. 

XT  Temporary  array  of  Interpolated  X values  for  crossings  at  level  H 
In  the  horizontal  direction. 

Y Table  of  meridional  grid  points  In  degrees. 

YF  Meridional  value  of  first  point  of  a string. 

YMIN  Minimum  value  of  Y after  scaling  to  10  Inches. 

YSYMB  Vertical  coordinate  of  leftmost  comer  of  plotted  character. 

YT  Temporary  array  of  Interpolated  9 values  for  crossings  at  level  H In 
the  vertical  direction. 

Z Table  of  functional  values  at  grid  points  X,  Y. 

^ ^max“^ln 

ZF  Zd/2*NCL  - AZ 

ZMAX  Maximum  value  In  Z table. 

ZMIM  Minimum  value  In  Z table. 

SUBROUTINE  DPFT 

Returns  the  2-dlmenslonal  spectrum  of  the  data  field  FLD  (containing 
the  variable  specified  in  DISPLY  at  a fixed  level) . A restriction  is  that 
all  latitudes  should  contain  the  same  number  of  grid  points  (NBJX). 
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Common  Blocks 


/DPC/ 

Variable  Description 


C Temporary  variable  holding  odd  components  of  north-south  sweep. 
CZ  Data  array  or  spectral  array  returned  by  FFT. 

FACl 

FAC2 


FAC3 


Normalizing  factors 


FAC4 


LA  NBJX/2  +1+1. 

IX  FFT  code. 

NT  2JB  set  in  DISPLY. 

N2  NBJX/2  set  in  DISPLY. 

N2P  N2  + 1. 


S Temporary  variable  holding  even  components  of  N-S  sweep. 


SV  Zj(2). 


The  logic  of  the  loops  is: 


J - 1,JB 
I - l,n 

Z±  ■ FyCi)  Fill  CZ  with  data  from  band  J 
6 

Z * FFT(Z)  Replace  data  in  CZ  with  E-W  spectral  components 
I - 2,2JB 

Fj(i)  - 4/n  Z odd 
Fj(i)  ■ -4/n  Z even 
7 

Fj(tt+2)  - (l/n)Z2 
Fj(l)  - (l/n)Zi 

Fj(n+2)  - Fj(n+2)  - Fj((n/2)+2)  - 0 
1 


125 


•I  - l,(n/2)+l 
•J  - 1,JB 


EJ  ■ 1,jB 

Zj  • Fj(j)  Fill  cz  with  E-W  components  i. 

Zl+JB  ■ Augment  data  using  Z(n/2)+j  = 2(n/2)-j;  which 

3^  for  real  data  -<«  Z(ii/2)+j  - 

Z “ FFT(Z)  Replace  data  in  cz  with  N-S  sweep  of  E-W  components, 


Z “ FFT(Z)  Replace  data  in  cz 

CJ  ■ 2,JB 

Fi(j)  - A/n  Z odd 

^I(JBl+j)  ‘ ^ 

4 

5 ^I(JBl)  ■’  5 Fi(JB2)  ■ ® 


- -4/n  Z even 


SUBROUTIIJE  FFT(CZ.M.KIV.IRE.IX) 

Returns  the  spectrum  of  the  complex  data  array  CZ-.  The  transformed 
data  replaces  the  original  data  in  CZ.  If  the  original  data  is  a real 
sequence  {xj}  j ■ 0,  1,  . . . ,n-l  we  must  have  n ■ 2*^,  with  k ^ 2.  Then 

(n/2)-l 

Xj  • 9.^/2  + Z [ait  co8(2irjk/n)  + bi^  sin(2irjk/nQ 
k-1 

f*n/2)^2  cos(iTj)  J • 0,  1,.  . n-1 

n-1 

with  ai,  • 2/n  Z k.  coa(2itjk/n)  k ■ 0,  1,.  . n/2 

j-0  ^ 

n-1 

bt,  “ 2/n  Z Xa  sin(2iTjk/n)  k ■ 1,  2,.  . n/2  - 1. 

J-0  ^ 

In  the  more  general  complex  case 

. - L ZTTlk]  . _ „ , 


*J  - Jq“x«"PP^J  ' 

where  is  determined  by  the  inverse  transformation 


n-1  r , 

Oi.  ■ 1/n  Z z.  exp  ii  - 
j-0  ^ ^ 

ao/2. 


lC*0f  lye  • ey  Q"!  e 


■ni.o  • ao/2.  3(^/2)  - a(„/2) 
^ 2 


k-1,  2,.  . .,n/2-l. 


•v  “ 

«ti-k  • 2 


For  the  forward  transform  (INV  » -1),  and  real  data  (IRE  ■ -1) 


where  t. 


(zo»  ^n/2^ 


^(n/2)-]! 
*(n/2)+l 

I-  *n-l  -i 


Returned 


L.^  - exp  Tl=-'^"  I - , and  for  1 < n < n-1  is  complex.  Since  Zg  and 
jit  ^ n j j — 

*n/2  real,  we  need  only  return  the  first  (n/2)  + 1 coefficients  on  top 
of  the  original  data,  by  using  the  following  method.  Zq  and  2^/2  occupy 
the  first  complex  number  (2  words).  Then  z^  thru  Z(n/2)-i  occupy  the  next 
n-2  words,  replacing  the  original  real  n data  words,  with  the  transformed 
values.  Then  since 


n-1 


Z „ 2iT(n/2  + r)k|  _ 

‘(n/2)+r  “k  n J 

- I exp  [i  - Z(n/2)-r 


1,  2,.  . .,  n/2  - 1 


k-0 

since  for  real  data  > o^. 
Then,  returning  Zg,  zj^,. 


» *(n/2)-l*  *n/2 


we  can  recover  Z(n/2)+i *n-l  *(n/2)+r  “ *(n/2)-r 

For  the  more  general  complex  case  • Sg.  For  the  forward  trans- 
form with  real  data  then,  it  is  enough  to  sequence  the  even  indexed'  data 
as  the  real  part,  and  the  odd  indexed  data  as  the  Imaginary  part  of  CZ. 


For  the  inverse  transform  (INV  ■ +1)  and  complex  data  (IRE  ■ +1) 
n-1 

k • 0,  1, . . . , n-1 
127- 


“k  ■ » ^ 


“1 


an-Li 

, 1 _*-► 
so  that  o ■ n T z. 

Then  the  values  returned  by  the  forward  transform  should  be  rearranged 

- ^ 

and  augmented  using  normalization  by  — is  not  done  in  FFT 


”*0 

n 

« 

T* 

• 

• 

Zn-l- 

itself. 


Buffer  out 
lOT  band 
by  band 


I 


In  the  following,  JUMP  locates  the  beginning  of  data  for  the  current 
band  (each  band  holding  384  data  words),  ILL  or  IGG  locates  Che  beginning 
of  data  for  a cell  (each  cell  holding  12  data  words),  and  IL  or  IG  locates 
the  current  variable  within  Che  cell.  HLD  is  dimensioned  128  to  hold  a 
full  filtered  baH^  of  4 variables  (128  > 4 x 32).  These  values  are  held, 
rather  cj|pn  j^edlately  put  on  Cop  of  the  old  values  in  BLK,  since  as  we 
go  around  Che  band  we  will  need  Che  old  values  of  I«  1,  2,  3,  4 for  Che 
9 point  filtered  value  of  I =■  32,  for  example. 

The  E-W  filter  is  straight  forward.  The  N-S  filter  requires  some 
ingenuity  to  conserve  core.  The  strategy  decided  upon  was  to  calculate 
the  9 point  filter  for  the  16  bands  of  fixed  zonal  index  I in  the  following 
pattern: 


1 point 


3 point  5 point  7 point 

4 


15  •- 

16  e- 


1 

2 

M . 

3 

1 

3 

6 

2 

5 

9 

4 

8 

13 

7 

12 

17 

11 

16 

ii 

15 

20 

25 

19 

24 

29 

23 

28 

.33 

27 

32 

_J7 

31 

36 

41 

35 

40 

45 

}9 

44 

49 

43 

48 

53 

47 

52 

57 

51 

56 

60 

55 

_ _ 

59 

. _ 

62 

50  y 

58  2 

63  7 
647 


Data  la  oavcd  in  an  array  indicating  the  various  stages  of  its  processing 

3 pt.  , 5 pt.,  7 pt.,  and,  when  it  reaches  the  9 pt.  stage  on  the  right  of 

Che  diagram,  it  replaces  the  original  data  on  the  same  level  in  the  diagram. 
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J-MMax 


Call 

DELSQ(JTJMP) 


I 


When  the  circled  number  (H)  reaches  64,  the  process  begins  again  with 


data  from  the  next  zonal  Index  I + 1.  JUMP  locates  cells  of  that  zonal 


Index  within  each  band,  with  excessive  data  values  separated  by  384  words. 


9 COS04.U  C0Sd44;< 


where  the  coefficients  CJHM 


cos64_l 


and  CJHP 


coa9.j^j^ 


have  been 


calculated  In  GEX,  and  stored  In  the  array  CF. 


DELSQ  returns  6^  of  the  variables  In  position  M-  successive  calls  to 


DELSQ  yield: 


for  N - 1 


5^gj  where  gj  Is  held  In 


for  N = 2 


2 2 2 

6^(5  gj)  where  £ gj  Is  held  In  the- 3 pt.  arrays 


for  N - 3 


for  IJ  ■ 4 


2 4 4 

£ (£  gj)  where  £ gj  Is  held  In  the  5 pt.  arrays 

gpj  ■ Q - gj  where  £®gj  Is  held  In  the 


7 pt.  array,  and  gp^  replaces  the  values  of  gj  In 


the  array  BLK. 


For  example,  the  first  10  calls  to  DELSQ  would  produce: 


t 


Values  out  of  array  bounds,  as  for  example  BLX(JSPl)  for  J ■■  16,  are  multiplied  by  Che  coefficient  zerol 


tFONfCHllOr oG  ,T3.  .’■WC  STtP  KO-i  nc  FILTER.  2j17  LETNfliiO 

ATMCH.GuDI-^,  C.iY2q7'i,IJ=,iETTA'icR,'<<  = L . 

ATTACH,'’‘aai-» , T.^,?;3r,A2R7'4,:D=M£T'^a  '£R,  MR=i. 

ATTACH,T4  = f3, 5R*?E*,GX2=}7V,  TOs-’iTTa  .l£<,  mR=1. 

RE3'JE^T»TAPE1,»3- . 

UPOATE  (3 

ETN,OPT  = t,i  = r;CMhILI,  SL. 

FTN,OPT=2,Su. 

pAusE.'GPi'^arcR  ~np  on 

LDSET(3?=.Si.T  = 7cR')) 



CATALOG.TIoET  TiMESXag-’t.,  lJ  = M£TTAtJER,H-<=x,RP=9^R. 

♦IDENT  -^SET 
♦0  bGEX.  3 

IT  A PEft  = '*  , T33£A  = : ,Ta3El  j = 0,rAP-:i...  = « jTAf’ilSsO  ,TAP£12  = J ) 

•I  i»r,Fx.6 

COMMCN/^OEF/OFCSS) 

CO,'i>ION/FTLT/J\/aK(4i  .NUM'’ 

■■■  C0'1M0N/0V?"3/0;ig'3,CJ 
OT-IEN'^I  ON  (wE'’) 

OiMErJSTCN  Tl^Cj) 

eguivalenc-;  (Clk(1)  , (G.i^  (c.)  ,n''ave) 

♦'I  liGEX.ll 

OA'^A  inFS2/l-'/ 

•I  GGtX.x6 

REWTNG 

*1"GGFX.2= 

PEAKTNPL'T.IEJ)  IFX  , IFlL.NJMg,  (JV3R(I)  ,r  = l,NUMR) 

■■■'ISG  FOPMSTCi  I5/,i5) 

MRIfFdP-lT,  1>1)  try,  tfIL.M'H3,  ( J'MR(T)  , In  j'l'JNP) 

' I5l  format  (xHu, 5.  <,*FTLT£r<  o 3 ? A «£  i £ Ro  * / //*  OO'^E  FOR  LAjT  TTM^  STEP  IF^ 
IP^'+XjxS/*  FREOUENPy  OS’  FTLiiv  x I =*  yX  , I 5/ ♦ NJMP=-rc  Or  VARiapLES  T"* 
2 filter  NUM3:*Iy/*  0^0^  '/ A R JC  GL  £ 3*  / * U = l*/*  \l=Z*/*  T = 3*/*  R=4* 

_ 3/*  osr=r*//4i3) 

ALOHh  = . £ *(P3vJn(2  )-PAN0(3  ) ) 

00  IfO  J=2,i’ 

J0=2*J-c  . JE=JO+i 

JPlsJri 

BANOj  = RArtC(  J)  33ANJP^  = 3A,4n(  j = : ) 

CP  (JO)  =C03(  OAMD  Ji  /CTriJ(j>  -C"  ( JE)  =C0S(3AN  Jr  l ) /CTH  J ( J) 

“ CONTIM'JF 

CJJ=-x.*COS(ALPHA) 

"'CP(1)=''FC».)  sj. 

0\/25c  = l./2"S. 

♦I  GGEX.  <i9 

N<V=K/*K''‘«=<'^-x 

*I  GGFX.3J  rO  "\ 

N'^tScSy  V U'  ^ 

NP  = ls,x3l6l3l61Ei61MPx6‘'  % ! ^ \\  L\ J i..'w  ^ 

IF  = jL46L6y5J  GOe  36j60oj<^'’  KM  iX\ 

msk=77occcl:.  otOcC'rofc.9 


r rO^H 

> r.  v^oi  ' 


irjj  r 

ii>^L  Lw'i  I 


♦I 

— BEST  miUZi  CC?Y 

_*^(i»;£X.a2 

N1.cx=tfx  ‘ . . 

14j  CONTI'lUt. 

■ '-TF  CIOTIM  S,  I- TL)  .EQ.  Q>  '3  TD  ■»? 

.is  (iN  + i.- ) /t.^ 

iSHFTsic  ♦MOn(^4,  IJ' 

IlsMXTFl  X(15<  ,IinF-^)  _ _ • 

i3=I*.4Mn..^s-  ■• 

lAf'xaYU  ) =ia\SAYf  m;  ♦T.i 
N=N>1 

i-X=ri30  T3  39 
37  CO‘'TIiUE 
IFX  = 1 

■ .is  (M+I  J ) /I, 

ISHFTsc,^  + (-5)  ♦MO''(N,Ij) 

Ils-iXis-i  X(iS<,  IS-IFT) 

I2=Il.ArO.C- 

I^o.^.MY  (^')  = IA'’RA''f i-it. 

N = M+1 

REWI.N3  lor' 

UisUiMlT  (lOT)  ^9UFFE^  IN(Tor,i)  ( :l  K ( 1) , CUX  ( 3 ) ) 

DO  5C1  Jsl,i.6 
UNsUMTT  (137) 

J2=oP4»J 

Jt=j4-3e3 

' ■'  9JFpr^  iNdDT,!)  (fl<  ( »1)  , 3.  >' ( J 2 ) ) 

50j  ''ONTIiUt 

"(UsUIJIT  (lOT) 

CAlL  FiLTc-W 
CALL  FILTnS 
.^£WT^)3  TOT 

" UisHNI-^  (Tcr)  ^9  iF^ER  -'UTdOTjL)  (C  L<  ( 1)  , CL<  ( 3 ) ) 

03  >C1  JsijiF 
■■  ■■J2=5  6^*J 

J*=J2-o83 

U.vIsUNTT  (TOT) 

0Jr'?E<  OUTdOFjl)  (9i.K'(  Ji.)  Ji)  ) 

“ 501  CQNTiiUE 

UNsUNTT  (TOT) 
no  •»  i = i ,3 

3 CALL  =JUFT(rjT,L,'’''s,'J,0L<,7)  ?CALl  CRJFT  (IOI  , 4, 9u9,M,  3LK  , 3' 

REWIH3  loT 

IF  (Ml'^X.EC.  1)  “»r»  "^3  3'^ 

REW^MO  TNi. 

U.isUNIT  (IN*)  *9'«F'^ER  I^'  (^NliL)  ( CK  ( 1)  , CLK  (3  ) ) 

03  Si'”  J=i»t' 
j2=384» J 

JlsJ2-seo  ... 

UNs'iNTT 

9UFFF.R  IN  (!»•!, 11  (^L  < ( »1}  , "LX  ( J2  ) ) 
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555  CON-^IMUfe 

'J'4='IMIT  (IMl) 

CALL  PILTE'J 

'CALL  -i(  TN5 
•tiWlNn  INx 

IJNSUMIT  (T'U)  '.3  te-f  E'»  QMTdNIi,!'  (rL<  ( *)  , CL<  (3  ) ) 

DC  S'*  2 J=l,*6 
J2=3P4* J 

_Jl  = J2-3<'5 

UN  = UHIT  (I'Jti 

a'JFe’F.<  OU'^d'Jl,!)  I aL<(Jl)  , 3.  <f  J2)  ) 

■ 5u2  CONTlNlUt 

UN=UNIT (IHL» 

0j"4"I*l,3 

4 CALL  ’3l'FT(Txii.  ,IJ,plK,3)  SCALL  C3JFT  ( I N1 , 4 * iJ » U,  TLK  , 3) 

^iWIHC  IN^ 

36  COMTiMME 

t?  (H00(NTS,  .Nl.  r>  ro  T3  35 

i- (NTS.Mt.MTi)  jO  TO  4w 
(INI) 

CALL  COh.P  (I'Jl,  i neS2  ) 

"E'luFiLt  I"!  52 

N3ESt  = ‘irE52*-L  ••  _ . _ 

REWIND  INI 
4^  CDNTTMUF 

UNsNNIT (Tn-) 

Call  CONO{iD1  ,iD£^^) 

" "P.'iOFT^c-  THESi 

• NDtSi=NrES6<-l. 

REWIND  TOT 
♦0  GGEX.4’f5r,Ey.>5 

39  I- (NT3.lt.  ITS)  GO  TO  3L 

W,<ITir(£p;^’-,  i IAf^RAV(<j  ,<  = l,A) 

103  F0P,NAT(1H.,12A1.  ) 

WRITc(  ) iO-5,N'>E5?,''..pA,I‘^y''?£'^I.'tO  IHEa  iREWINDT0tS2'ST'^P 

' ion  pOkIAT  (♦  OIDTAL  <ri;.£5  D^iNF.RATED  *,2I5/»r  ''ODES  *215) 

•COMPILE  GG£X»GTTS.F,G3TTc-Lr 

»CDHPILr‘GF'<ri,GT’,GPSY*,  j'5o..,oPSTi.  oycLi 

•COMPILE  'AiOTATE,  n d,  v^ENTIS,:,^-  ,G3'J- 
“ •COMPILE  GSTiSF,GT£3. 

RJ’RuJTINE  -£LT=.n 

COMMO  J/FMF/9L<(  6«:23) 

CDMMC'l/FILr/  )'/hO(4)  .NMMO 
CDMMOM/CVER'i/Oy.iEo.  ''J 

OIME‘»SION  ML'^(129) 

'V  ■ J = 0 5J'JMFi-3'^'^ 

1 CONIICUL 
'HaO 

« JsJ+i 

I-(J.ED.l')  iD  TO  5 
JiJMEsJUMF^-OCL 

T-l'isES'IMSs’CNf  12  = 3 j^lMl  = 31  I = 32«  T=>i=l If-c  = 2 ‘IP  3=  3“.  TP  A = 4 


BESr  AVAIMBIE  COPY 


2 CONTINUE 

IN4=I  13f  Iri3  = r N?$I'12  = TN15  T>«1=I^I  = Ioi-IPi  = T0  2SIP2  = I03SI'03=I04 
_ i04  = '^  = 4+l 

Ir  (IFf.EO.33)  104=1 

^[LLf=JUfPf(  rN4-i)  *12^^U5  = JLMo+(IM3-i)  *12 

rLL2  = JJ''F*(r-t2-i)  *1  2-ILLl  = J'JM-»f(  iMi-i)  f^ablVVs  JU^nf  (I-l)  *12 

IOGl=  JUMP  + { [•'l-D  ♦12’‘Xo':?=  JU1P>  (iP.;-*)  ♦ i 2S  = J'MH*  ( lo-*-!)  *12 

I0G4r  JUMPf(  ro4-l)  *.2 
_ _IA  = a 

"3"  CONTINUE  " - 

IF (IA.EO.NJ  NO)  GO  I 0 4 

IA  = I/'fl 
IVAo=  JVAK  (lO 

lL4=lLL4fIVAPSTL3=TLL‘»fIVA^;£-?  = TLL2+I\/A'<‘‘ILl  = H.i.xf  TNAo 
I7=IV/f  IVAO.  rGl  = IG31fTvsp,X3  2=i--,(i2fI'/AKiIG3=i.GG3>I7rtR 
IGf=IGG4f i; AR  ” 

H=IHfl 

Hl'ndri)  = 07203*  (!_  ac,,  •^LK(IV)-l3LK(  ILfX-OLKC  IG4)  ) + 8.*  (OLXdLO)* 

19L<(IG3))-?«.*(Oi.v’(TL^jfiLi<(T".2))+G6.*(iL<(ILl)*'’i_<(lG*))) 

'“GO  TO  3 

4__IO(I.£0.32)  go  to  6 

CO  TO  2 ■ ■ . “ ‘ 

5 CONTI  (ME 

RETURN  “ 

6 .CON’TvjUE 
H=o 

_ II=JUNo-12 
CD  7 1 = 1, o2 
_II=I1+12 
IA  = 0 

continue 

IF  (lA.EC.iNJ  N'l)  CO  TO  7 

TA  = I''  + 1 
TVA-Js  J7AF  (lA) 

LLsIT+iVAR 

"iHslPfi 
9.K{LL)  = ULD{£rt) 

GO  TO  9 " 

7 CONTINUE 
GO  TO  1 

ENO 

SJBRf'J’i  TNL  FILTNS 

CJN'lON/rSO'Cjn  JMl,  JNEtOPl 

COMMON/ FMp /PL  <(65<;3) 

COMMON/ RO-T £/ UJH T (=)  .mcP’"  ( J)  , J7'»t  ( o)  .v^PK  o)  ,<?pT(3)  ,'f7PT(3)  , 

“ 1T3 PT  (3)  , T^artl)  ,T7=T  , atpr  ( I)  ,5  CRT  (5)  , i/PT  r31  ,PST3'’T(3)  jPST^nTd 

2),PoT7ot(3) 
k,OMMON/^Od* /OF  ( 

CJMMCJ/Fi!  T/J'/AR(4)  ,NU**P 

irOMMON/CV2  5 j/072  3C,CJ 

xauaJUMPs-^:  . * - 
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BEST  AVAILABLE  COPY 


V'^ONTlMUt 

JU'lPs  JjMPtL2 
m=0-‘M-IAX  = . 

T'CONTI'JUF 

M'iAX=MMflY4-i. 

^sHMAX 

_ _Ir  (J.'‘,T.16)  Sin" 

■ ^=1'  ■ 

3_C0NTI?'IUt 

'•l-sM+l  “ . 

J^l  = J-ii  J‘^2  = J-2’;  JP1  = J+i  I 

■— JTHM=2*  J-1*  jth?  = jp-im+;.  I 

CJHM  = 3F  ( S:JHP=^F  (JTH3)  ^ 

‘CALL  D£LSO(JU'4ci 
I-(M.EQ.o,)  *0  TT  4 

J=J-x  

NsN+l 

J JPsMCC (N,5) *J 
TF  (JU=.Fn.,  ) G'>  TO  ^ 

■■'■"GO  TO  3 

4 I-(I.L‘^.32)  GO  Tn  . 

■■■‘2FTUP.M 

5 CON'^T'JUE 
■ J=l5 

N=M1AX-15 

GO  TO  0 
£*^0 

SOBROOTIMt-  DiLSO(J'lMP) 

C0M'1CUCS'l/0J«4,''JHO,  rn,  Ji2»  JPl 
" COMMO'l/t'U‘^>'»wK(  6=5tS) 


COM'lLN/kOUT£/'II.  PT  (oJ,'J-3T(3»,j7PT(3>,'npT(3),\/5=T(3)  , 7 7PT  (3)  , 

1T3PT(3)  ,TF5r(f),»70T(3),13  = r(3)  ,R5=>T(3)  , P7  PT  ( 3 ),  PST  3’^T  ( 3 ) ,PST5PT  (3 

3\  Q^T*»n*r  / 


2) ,PST70Tfo» 

■ OOMMON/FILT' J'/A^s  r4)  ,M"mq 
C0mv|0M/CV25o/O  Vt  F6,  r J 
■■  OIMtNilON  A3-.AV(L5) 
cOUIVALtNCi  ( Ar,<AV(l)  ,'J3»r  (1)  ) 
lAsQ 

IF { (N-1 ) • (N-4) ' 

: CON^lMttE 

_J1J  VsMOO  (J 

MJVisiuC  (JM*. , 7) 

.lJV2=‘^OOt  J*t’,3' 

N-12TTa5*  (M-2> 

IF  (IA.EC.MJ'11)  => 
lAsIAf^ 

IVA,<=J\/An.(l  A) 

rNC=s*  a n 

_L'1UVA  ?sI 
J3sL0CVA,<+‘»071 
J-<=JF>3 


:iUr<.N 


\ f 


Ml  *0^  r 

AiLrAiiLL 


COrY 


JSPl=LOCVfl^*-'^OV  . 

J3Ml=LOC'»a<<-'^OV/.: 

AROAV(JR)  =:  J-iM*ACRA  Y ( JSMl)  f : j»aP'(AY  f ,J3)  +CJHP*  » JS«1) 

';'D  TO  .i 
~IJ  ■ CONTINMF 

M0VBs384»JMi«- JUMP 

flovis^orc 

I"(.'4.£a.t)  "3  TO  5»!j 
■■"■3  rONTiNUt 

IFlIA.EC.N'JilP'  PiTU^M 
* lAsIA+l 

IJA?3;JVAa(IA) 

INC=S*(  IV'^v-lWi 

JSalVAi^+MOi/’ 

J3'Mi=JS-3i'^5IPU.EM.l)  J3Mi  = JS 

J3BlrJS^3S4 

LOCRalMC 

ARRAY  (JR)  =:  J-tM*PL<(  J3M1)  4-0  J* '’L  K ( JS)  ♦CJH  P’BLK  ( JSPl ) 

GJ  TO  3 

— <7~C0NTI,<JUF  ■■  ■" 

MOV  = MOO(J,i) 

■ MOVca-'On  (j<?, 

IF  (lA.tC.'MMd)  OcTU^N 

TAalA+i 

IJ  ARaJVARd  A) 

IMC=«*( TVA3-t) tx 

JRsrVA-<+HO'/» 

t.DCVAR=IHC<-5 

JSsLCCV  AR+HO'/l 

JROiscOCV  A^if  10V 

J>Mi  = L0CVA=  ♦■I  J V’ 

r6L<(J^)  -0V253*  (wJHM*  ARRAY  ( JSMi)  ♦rj*A(^P.AV(JF)  + OJHP* 
__  *A^RAY{  JSF*)  ) 

GJ  tq  2C 
E'lO 

R=4  LF  TtST  FIlTcX 

S3  ARAM  FTI  = 36-.»  >iTi  = 2 , T^cSFa^  , *•  TT  = 72  0 . , RTL=*»0 . 3 


For  the  codes  below  1 < N < NBLK 


2 <.  J <.  JBl 
1 < L < 4 


SUBROUTINE  CODES 


INIT 

(15.1) 

(15,2) 

(16,0) 

(16,1) 

OBAL 

(20 

.N) 

REPORT 

(30,1) 

(31, J) 

(32.J) 

(33,J) 

(34,0) 

STIS4» 

(95,0) 

(+60,0) 

(±61,0) 

(+62,0) 

(±63,0) 

(±64,0) 

TIS<^ 

(95.1) 

(±55,5) 

(±56, J) 

GEX 

(60,0) 

(60,1) 

COMP 

(+1.0) 

(+1.1) 

(±3.J) 

(±4.N) 

ENTRY  EXPAND 

(±2,1) 

(+5,N) 

(±6,J) 

(±7,L) 

ENTIS 

(99,0) 

(99, L) 

BDUMP 

(40,0) 

(40,+l) 

(40,2) 

(40,3) 

DISPLAY 

(0,0)  (0,N) 

SKFL 

(1,0)  (1,N) 

SKFL 


